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Abstract 

The basics of the relativistic astrophysics including the celestial 
mechanics in weak field, black holes and cosmological models are illus- 
trated and analyzed by means of Maple 6 

Application Areas/Subjects: Science, Astrophysics, General Relativity, Ten- 
sor Analysis, Differential geometry, Differential equations 
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1 Introduction 



A rapid progress of the observational astrophysics, which resulted from the 
active use of orbital telescopes, essentially intensifies the astrophysical re- 
searches at the last decade and allows to choose the more definite direc- 
tions of further investigations. At the same time, the development of high- 
performance computers advances in the numerical astrophysics and cosmol- 
ogy. Against a background of these achievements, there is the renascence of 
analytical and semi-analytical approaches, which is induced by new genera- 
tion of high-efficient computer algebra systems. 

Here we present the pedagogical introduction to relativistic astrophysics 
and cosmology, which is based on computational and graphical resources 
of Maple 6. The pedagogical aims define the use only standard functions 
despite the fact that there are the powerful General Relativity (GR) oriented 
extensions like GRTensor |lj . The knowledge of basics of GR and differential 
geometry is supposed. It should be noted, that our choice of metric signature 
(+2) governs the definitions of Lagrangians and energy-momentum tensors. 

The computations in this worksheet take about of 6 min of CPU time (PIII- 
500) and 9 Mb of memory. 

2 Relativistic celestial mechanics in weak gravita- 
tional field 

2.1 Introduction 

The first results in the GR-theory were obtained without exact knowledge 
of the field equations (Einstein's equations for space-time geometry). The 
leading idea was the equivalence principle based on the equality of iner- 
tial and gravitational masses. We will demonstrate here that the natural 
consequences of this principle are the Schwarzschild metric and the basic 
experimental effects of GR-theory in weak gravitational field, i. e. planet's 
orbit precession and light ray deflection (see Q). 

2.2 Schwarzschild metric 

Let us consider the centrally symmetric gravitational field, which is produced 
by mass M. The small cell falls along x-axis on the central mass. In the 
agreement with the equivalence principle, the uniformly accelerated motion 
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locally compensates the gravitational force hence there is no a gravitational 
field in the free falling system Koo- This results in the locally Lorenzian 
metric with linear element (c is the velocity of light): 

d s 2 = d Xoo 2 + d yoo 2 + d z^ 2 - c 2 d too 2 

The velocity v and radial coordinate r are measured in the spherical system 
K, which are connected with central mass. It is natural, the observer in 
this motionless system "feels" the gravitational field. Since the first sys- 
tem moves relatively second one there are the following relations between 
coordinates: 

^=7=1 (f* = V 



d too = yr 3 ^ dt 



d Uoo =rd 6 



d Zoo =t sin(0) d 4> 

The first and second relations are the Lorentzian length shortening and time 
slowing down in the moving system. As result, Koo from K looks as: 

(is 2 = (1- Z? 2 )^ 1 ) dr 2 +r 2 (d 6 2 + sin(#) 2 d (j) 2 ) - c 2 (1 - (i 2 )d t 2 

The sense of the additional terms in metric has to connect with the char- 
acteristics of gravitational field. What is the energy of Koo in K? If the 
mass of Koo is m, and mo is the rest mass, the sum of kinetic and potential 
energies is: 
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> restart: 

> with(plots) : 

> (m - m0)*c~2 - G*M*m/r=0;#energy conservation law\ 

> (we suppose that the Newtonian law of gravitation\ 

> is correct in the first approximation) , G is\ 

> the gravitational constant 

> °/./(m*c-2): 

> subs( m=mO/sqrt (l-beta~2) ,% ):# relativistic mass 

> expand (°/ ) ; 

> solve( 7o,sqrt(l-beta~2) ): 

> sqrt (l-beta~2) = expand (%) ; 

> l-beta~2 = 

> taylor ( ( 1-subs (op (2 , %) =alpha/r , rhs (%) ) ) "2 , alpha=0 , 2) ; \ 

> #alpha=G*M/c~2 , we use the first-order approximation \ 

> on alpha 



(m — mO) c 



2 GMm 



l- x /l-(3 2 -^ = 



GM 

c 2 r 



1-/3 2 = 1 



GM 

c 2 r 



l-/3 2 = l-2-a + 0(a 2 ) 
r 



The last results in: 

2 _ dr 2 , J2 i j eft , -/p\2 j jjl \ „2 / 1 _ ga 



S 



+ r 2 (d 9 2 + sm(6) 2 d (j) 1 ) - c 2 (1 - ?f )d t 2 



This linear element describes the so-called Schwarzschild metric . In the 
first-order approximation: 

> taylor ( 

> d(r)~2/(l-2*alpha/r)+r~2*(d(theta)~2+sin(theta) ~2*\ 

> d(phi)~2)-c~2*(l-2*alpha/r)*d(t)~2, alpha=0, 2 ): 

> convert ( °/„ , polynom ) : 

> metric := d(s)~2 = collect( °/„ , {d(r)~2, d(t)~2} ); 

metric := d(s) 2 = 

2 + 2c^a d 2 + , 1+ 2a, d , x 2 + r 2 ™ 2 + gin ^)2 d /0N2) 

r r 
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2.3 Equations of motion 



Let us consider a motion of the small unit mass in Newtonian gravitational 
potential $ = - . Lagrangian describing the motion in centrally sym- 
metric field is: 

> L :=\ 

> (diff (r(t) ,t)~2 + r(t)"2*diff (theta(t) ,t)~2)/2\ 

> + G*M/r(t) ; 

1 r 9 , XN 9 1 , xo , d „. ,.n GM 



The transition to Schwarzschild metric transforms the time and radial differ- 
entials of coordinates (see, for example, [||, Q): dr > (1 + ^ )dr and dt -> 
(1 + f-)^ 1 ** dt (it should be noted, that the replacement will be performed 
for differentials, not coordinates, and we use the weak field approximation 
for square-rooting): 

> L_n := collect(\ 

> expand (\ 

> subs(\ 

>{diff (r(t) ,t)=gamma(r(t))~2*diff (r(t) ,t) ,\ 

> diff (theta(t) ,t)=gamma(r (t) ) *dif f (theta(t) ,t)}, L ) ),\ 
>{diff (r(t) ,t)~2, diff (theta(t) ,t)~2}) ;# modified\ 

> Lagrangian, gamma (r) = 1+alpha/r (t) 

L-n := \ 7 (r(i)) 4 (| r(t)f + \ At? iW)? (| 0(t)f + ™ 



Next step for the obtaining of the equations of motion from Lagrangian is 
the calculation of the force ^ L(x, y) and the momentum L(x, y) , where 

y=m x: 
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> el := Diff (Lagrangian(r, Diff(r,t)), r) = \ 

> diff (subs(r(t) = r,L_n) ,r) ;#f irst component of force 

> e2 := Diff (Lagrangian(r, Diff(r,t)), Diff(r,t)) = \ 

> subs(x=diff (r(t) ,t) , diff (subs (diff (r(t) ,t)=x, L_n),\ 

> x));#first component of momentum 

> e3 := Dif f (Lagrangian(theta, Dif f (theta, t) ) , theta) 

> dif f (subs (theta(t) = theta, L_n) , theta) ; #second\ 

> component of force 

> e4 := Diff (Lagrangian(theta, Dif f (theta, t) ), \ 

> Diff (theta, t)) = subs(y=diff (theta(t) ,t) ,\ 

> diff (subs(diff (theta(t) ,t)=y, L_n) , y));#second\ 

> component of momentum 



9_ t — n „^ n „r„ A 
' at 



r 



el := jfc Lagrangian(r 
2 7W 3 (| r) 2 (£ 7 (r)) + r 7 (r) 2 (| 9{t)f + r 2 7 (r) (f 6(t)f ( £ 7 (r)) 



e2 := Diff(Lagrangian(r, —r), — r) = 7 (r(i)) 4 (— r(i)) 



e5 := A Lagrangian(0, | 0) = r(t) 2 7 (r(t)) 2 (| (9) 61) 



:= DifT(Lagrangian(0, -(9), -0) = r(t) 2 7 (r(t)) 2 (-0(t)) 



The equations of motion are the so-called Euler-Lagrange equations 

L(x, y) — L(s, y) = (in fact, these equations are the second New- 
ton's law and result from the law of least action). Now let us write the 
equations of motion in angular coordinates. Since e3 = due to an equality 
to zero of the mixed derivative, we have from the equation of motion 
0^ L(r, 9, yi, y 2 ) - ^ L(r, 6,y 1 ,y 2 ) = ( y x = §- t r , y 2 = §- t ) in the 
form: 

> Eu_Lagr_l := Diff (rhs(e4) ,t) = 0; 

Eu_Lagr_l := | r(t) 2 7 (r(t)) 2 (J^(t)) = 



Hence 
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m 2 im?ii t e{t)) = H= const 



> sol_l := Diff (theta(t) ,t) = solve(\ 

> gamma~2*diff (theta(t) ,t) /u(theta) "2 = H,\ 

> dif f (theta(t) ,t) );#u=l/r is the new variable 



sol A :=_ CW = _^L 



The introduced replacement u= j. leads to the next relations: 

> Diff(r(t), t) = diff (l/u(t) ,t) ; 

> Diff(r(t), t) = diff (l/u(theta) ,theta)*\ 

> dif f (theta(t) ,t) ; # change of variables 

> sol_2 := Diff(r(t), t) =\ 

> subs (diff (theta(t) ,t) = rhs(sol_l),\ 

> rhs(°/ ));# the result of substitution of\ 

> above obtained Euler-Lagrange equation 



» r(t) = _i^) 



(£u(0)) 



dt l{t) ~ u(0) 2 



2 



dt 7 



The last result will be used for the manipulations with second Euler-Lagrange 
equation L(r, 9, yi, y 2 ) - ^ L(r, 9, yi, y 2 ) = 0. We have for the right- 
hand side of e2: 

> Diff( subs( {dif f (r (t) ,t)=rhs (%) ,gamma(r (t) )=gamma} , \ 

> rhs(e2)) , t) ; 

|(-7 2 (|u W )tf) 



and can rewrite this expression: 
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> -2*gamma*H*diff (1+alpha/r (t) ,t)*diff (u(theta) ,theta) -\ 

> H*gamma~2*diff (u(theta) ,theta$2) *dif f (theta(t) ,t) ;#from\ 

> the previous expression, definition of gamma and\ 

> definition of derivative of product 

> first_term := subs({r(t)=l/u(theta) ,\ 

> diff (theta(t) ,t)=rhs(sol_l)},\ 

> subs(diff (r(t) ,t) = rhs(sol_2) //,)) ;# this is \ 

> a first term in second Euler-Lagrange equation 

a ^°(^»(*-W - g7 . ( ^ nW)( | W) 



first Jerm != - 2 *3*mfeW - H> < * »«,» u(Sf 

7 Ou 



el results in: 

> 2*gamma(r) ~3*diff (r(t) ,t)~2*diff (gamma(r) ,r) + \ 

> r*gamma(r) ~2*diff (theta(t) ,t)~2 + \ 

> r~2*gamma(r)*diff (theta(t) ,t)~2*diff (gamma(r) ,r) - \ 

> G*M/(r~2) ;# This is el 

> subs(\ 

>{diff (r(t) ,t) = rhs(sol_2), dif f (gamma(r) ,r) =\ 

> diff ( 1+alpha/r , r ) , \ 

> diff (theta(t) ,t) = rhs (sol_l) } , \ 

> °/ ):# we used the expressions for dif f (r (t) ,t) , gamma(r)\ 

> and the first equation of motion 

> subs (gamma(r) = gamma, °/ ) : 

> second_term := subs(r = l/u(theta) , °/ ) ; 

2 7 (r) 3 (| r(t)) 2 (& 7 (r)) + r 7 (r) 2 (| 6(t)f + r 2 7 (r) (| 6{t)f (& 7 (r)) 
GM 



second_term :- 



_ 2 ff 2 au(g) 2 (^u(g)) 2 + u(9fHi _ ff 2 u(g) 4 q _ QM<9f 



And finally: 



> Eu_Lagr_2 := expand( simplify (first_term -\ 

> second_term) /u(theta) ~2/H~2) ; 
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„ r ,d 2 u(0) n(6) 2 a GM 

Eu.Lagr.2 : = -(^ u(9)) - + + 



In the first-order approximation: 

> gamma~n = taylor ( (l+alpha*u) ~n, alpha=0,2); 

7 n = l + nua + 0(a 2 ) 



So 



> taylor ( subs(gamma = l+alpha*u(theta) ,Eu_Lagr_2) , \ 

> alpha=0 , 2 ) : 

> basic_equation := convert (%, polynom) = 0; 

q2 G M 

basic .equation := — u (®)) ~~ u (^) H — rpr + 3 u(6>) 2 a = 



2.4 Light ray deflection 

In the beginning the obtained equation will be used for the search of the 
light ray deflection in the vicinity of a star. The fundamental consideration 
has to be based on the condition of null geodesic line d s 2 =0 for light, but 
we simplify a problem and consider the trajectory of the particle moving 
from the infinity. In this case H= oo . 

> eq_def := subs(G*M/(H~2)=0,basic_equation) ; 

eq.def : = -(-^ u(0)) - u(0) + 3 u(9) 2 a = 



The free propagation ( a =0) results in: 

> subs (alpha=0 , lhs(eq_def) ) = 0; 

> sol := dsolve({%, u(0) = 1/R, D(u)(0) = 0}, u(theta)); 

> # theta is measured from perihelion, where r = R 
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;|l u( 0)) - U (0) = o 



cos(0) 

S0 Z := u (0) = — — 



That is r= co ^g) . The last expression corresponds to the straight ray 
passing through point 9 =0, r=R. To find the corrected solution in the 
gravitational field let substitute the obtained solution into eliminated term 
in eq_def: 

> -op(l,lhs(eq_def )) - op(2 , lhs (eq_def ) ) =\ 

> subs (u(theta)=rhs (sol) , op(3,lhs(eq_def ))) ; 

> sol := dsolve(\ 

>{°/„, u(0) = 1/R, D(u)(0) = 0}, u(theta)) ;#corrected\ 

> solution in the presence of field 



(-a + R) cos(0) 1 acos(2 6») o a 



The equation for asymptote lim^^ u(#) describes the observed ray. 
Then angle of ray deflection is ^ (symmetrical deflection before and after 
perihelion). Hence 

> simplify ( 2*subs ({theta=Pi/2 , alpha=G*M/c~2} , \ 

> rhs(sol))*R ) ; 

GM 



This is a correct expression for the light ray deflection in the gravitational 
field. For sun we have: 
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> subs({kappa=0.74e-28, M=2e33, R=6.96el0},\ 

> 4*kappa*M/R/4.848e-6) ; 

> # in ["] , where kappa = G/c~2 [cm/g] , \ 

> 4.848e-6 rad corresponds to 1" 



1.754485794 



The ray trajectory within Pluto's orbit distance is presented below: 

> K := (theta, alpha, R) -> 1 A 

> (- (alpha-R) *cos (theta) / (R~2) -l/2*alpha*cos (2*theta) A 

> (R~2)+3/2*alpha/(R~2)) : 

> S := theta -> K(theta, 2.125e-6, 1) :#deflected ray 

> SS := theta -> 1/cos (theta) : #ray without deflection 

> pi := polarplot(\ 

> [S , theta->theta , Pi/2 .. -Pi/2] , axes=boxed) : 

> p2 :=\ 

> polarplot( [SS,theta->theta,Pi/2. .-Pi/2] ,\ 

> axes=boxed, color=black) : 

> display(pl,p2,view=-10700. .10700 A 

> title =< deflection of light ray');\ 

> #distance of propagation corresponds to\ 

> 100 AU=1.5elO km, distance is normalyzed to Sun radius 



deflection of light ray 



10000 
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-5000 



-10000 







0.2 



0.4 



0.6 



0.8 
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2.5 Planet's perihelion motion 



Now we return to basic-equation. Without relativistic correction to the 
metric ( a =0) the solution is 

> sol :=\ 

> dsolve ({subs ({G*M/ (H~2) =k , alpha=0} , \ 

> basic_equation) ,u(0)=l/R,D(u) (0)=0},u (theta)) ; 



sol := u(8) = k 



(kR- l)cos(0) 
R 



This equation describes an elliptical orbit: 



u=k(l+e* cos(0) ), 



where e = ( ^ - 1) is the eccentricity. For Mercury k = 0.01, e 
0.2056, R = 83.3 

> K : = (theta, k, R) -> l/(k-(k*R-l)*cos(theta)/R) : 

> S := theta -> K(theta, .01, 83.3): 

> polarplot( [S,theta->theta,0. .2*Pi] ,\ 

> title =t 0rbit of Mercury'); 



Orbit of Mercury 
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The correction to this expression results from the substitution of ob- 
tained solution in basic-equation. 



> subs({G*M/(H~2)=k, 3*u(theta) ~2*alpha=\ 

> 3*alpha* (k* (l+e*cos (theta) ) ) ~2} , basic_equation) ; 

> dsolve({°/ ,u(0)=l/R,D(u) (0)=0} ,u(theta) ) ; 



-(^2 u(0)) - u(0) + k + 3 a k 2 (1 + ecos(#)) 2 = 



u(9) = 

3 1 

3ak 2 ecos{9) + 3sm(9)ak 2 e9 + -ak 2 e 2 - - a k 2 e 2 cos(2 9) + k + 3 a k 2 

(3 a k 2 e R + a k 2 e 2 R + k R + 3 a k 2 R - 1) cos(0) 

R 

Now it is possible to plot the corrected orbit (we choose the exaggerated 
parameters for demonstration of orbit rotation): 

> K : = (theta, k, R, alpha, e) -> 1 A 

> (3*alpha*k~2*e*cos (theta) + \ 

> 3*sin(theta)*alpha*k~2*e*theta + \ 

> 3/2*alpha*k~2*e~2 - l/2*alpha*k~2*e~2*cos (2*theta) +\ 

> k + 3*alpha*k~2 - (3*alpha*k"2*e*R + alpha*k"2*e"2*R +\ 

> k*R + 3*alpha*k~2*R - 1) *cos (theta) /R) : 

> S := theta -> K(theta, .42, 1.5, 0.01, 0.6): 

> polarplot( [S,theta->theta,0. .4.8*Pi] ,\ 

> title =< rotation of orbit'); 



rotation of orbit 
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Now we try to estimate the perihelion shift due to orbit rotation. Let's 
suppose that the searched solution of basic-equation differs from one in the 
plane space-time only due to ellipse rotation. The parameter describing this 
rotation is uj : u( 6 ) = k(l + e cos^ — ujO) ). 

> subs({G*M/(H~2)=k, u(theta)=\ 

> k*(l+e*cos (theta-omega*theta) )} , \ 

> basic_equation) : #substitution of approximate solution 

> simplify (%) : 

> lrisC/O: 

> collect (%, cos (-theta+omega*theta) ) : 

> coeff(%, cos (-theta+omega*theta) ) : #the coefficient\ 

> before this term gives algebraic equation for omega 

> subs(omega"2=0,7o) :#omega is the small value and \ 

> we don't consider the quadratic term 

> solve (7o = 0, omega): 

> subs({alpha = G*M/c~2, k = 1/R/ (1+e) } , 2*Pi*%) ;#result\ 

> is expressed through the minimal distance R between\ 

> planet and sun; 2*Pi corresponds to the transition to \ 

> circle frequency of rotation 

> subs (R=a* (1-e) ,%) ; #result expressed through larger\ 

> semiaxis a of an ellipse 

ttGM 
6 c 2 R(l + e) 

g ttGM 
6 c 2 a(l-e)(l + e) 

Hence for Mercury we have the perihelion shift during 100 years: 

> subs ({kappa=0 . 74e-28 ,M=2e33 , a=57 . 9ell , e=0 . 2056} , \ 

> 6*Pi*kappa*M/a/(l-e~2)*(100*365. 26/87. 97)/4.848e-6) :\ 

> #here we took into account the periods of Earth' s\ 

> and Mercury's rotations 

> evalf (%);#["] 

43.08704513 

Now we will consider the basic-equation in detail ||. The implicit solu- 
tions of this equation are: 
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> dsolve( subs (G*M/H~2=k, basic_equation) , u(theta)); 

MO) i 

/ = d_a-9- _C2 = 0, 

J V-Cl -_a 2 + 2k.a + 2_a 3 a 

MO) 

/ - — (La - - J'2 - 



MO) i 

yfjCl - -a 2 + 2k_a + 2_a 3 a 



Hence, these implicit solutions result from the following equation ( (3 is the 
constant depending on the initial conditions): 

> diff( u(theta) , theta)~2 =\ 

> 2*M*u(theta) ~3 - u(theta)"2 + 2*u(theta) *M/H~2 + beta;\ 

> # we use the units, where c=l, G=l (see\ 

> definition of the geometrical units in the next part) 

> f := rhs(%) : 

> f = 2*M*\ 

> ( u(theta) - ul )*( u(theta) - u2 )*( u(theta) - u3 ) ; 

(|u W ) 2 = 2Mu W 3-u(^ + ^§^ + /3 



2 M u(6f - u(d) 2 + 2U ^ M + (3 = 2 M (vl{9) - ul ) (u(0) - u2) (u(9) - u3) 



Here ul, u2, u3 are the roots of cubic equation, which describes the "poten- 
tial" defining the orbital motion: 

> fun := rhs(%) ; 

fun :=2M (u(0) - ul ) (u(0) - u2) {u(9) - u3) 



The dependence of this function on u leads to the different types of motion. 
The confined motion corresponds to an elliptical orbit 

> plot (subs ({u3=2, u2=l, ul=0.5, M=l/2} ,f un) ,u=0 . 4 . . 1 . 1 , \ 

> title= ' elliptical motion' , axes=boxed, view=0 . . . 08) ; 
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The next situation with u->0 (r-> oo ) corresponds to an infinite motion: 

> plot (subs ({u3=2, u2=l, ul=-0.1, M=l/2} ,f un) ,u=0 . . 1 . 1 , \ 

> title= 'hyperbolical motion' , axes=boxed, view=0 . . . 6) ; 



hyperbolical motion 




Now we return to the right-hand side of the modified basic equation. 
> f ; 
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In this expression, one can eliminate the second term by the substitution u( 

> collect ( subs (u(theta) =\ 

> (2/M)~(l/3)*y(theta) + 1/(6*M), f ) ,y (theta) ) ; 

3 



4y(#r+ 

t 2^/3) f_L)(l/3) 2 2(V3) (J_\(l/3) M \ 

6 M # 2 

V / 



This substitution reduced our equation to canonical form for the Weier- 
strass P function || : 

> g2 = simplif y(coef f (%, y(theta))); 

> g3 = -simplif y(coef f (%%, y (theta) , 0)); 

> diff( y (theta), theta) "2 =\ 

> 4*y(theta)~3 - g2*y(theta) - g3; 



92 = - 



1 2 ( 1 /3)(_L)(l/3) ( _ i? 2 + 12M 2 ) 



6 MH 2 

1 —H 2 + 54/3 M 2 H 2 + 18 M 2 
93 ~ "54 APH 2 

(-^y(O)) 2 = 4y(ef-g2y(0)-g3 

that results in: 

> y(theta) = WeierstrassP (theta, g2, g3) ; 

y (0) = WeierstrassP (9, g2, gS) 
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In the general case, the potential in the form of three-order polynomial 
produces the solution in the form of Jacobi sn-function: 

> Orbit := proc(f, x) 

> print ( 'Equation in the form: \i' (theta) "2 =\ 

> a [0] *u~3+a [1] *u~2+a [2] *u+a [3] ' ) : 

> degree (f ,x) : 

> if (°/. = 3) then 

> a[0] := coeff(f, x~3) :# coefficients of polynomial 

> a[l] := coeff (f , x~2) : 

> a [2] := coeff (f, x) : 

> a [3] := coeff (f, x, 0) : 

> sol := solve (f =0, x) : 

> print ('Roots of polynomial u[l] < u[2] < u[3]:'): 

> print (sol [1] , sol [2] , sol [3] ) : 

> solution := u[l] + (u[2] -u[l] )*JacobiSN(theta*sqrt(\ 

> 2*M*(u[3]-u[l] ) )/2 + delta ,\ 

> sqrt ( (u [2] -u [1] ) / (u [3] -u [1] ) ) ) "2 : 

> print( 'Result through Jacobi sn - function'): 

> print (u (theta) = solution): 

> else 

> print ('The polynomial degree is not 3') 

> fi 

> end: 

Equation in the form : 

u'(theta)"2 = a[0] * u" 3 + a[l] * u" 2 + a[2] * u + a[3] 
Roots of polynomial u[l] < u[2] < u[3] : 
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iCl/3) 



1 70 

6 MH 



1 -H 2 + 12M 2 e 
6 Mi? %l (1/3) + M' 



j_%l(V3) ^(-i? 2 + 12M 2 ) | 
12 Mi? + Mi?%l (1/3) + M + 



1 ^cl^ 3 ) 



6 Mi? 



+ 



(-i? 2 + 12 M 2 



V 



Mi?%l 



(1/3) 



1 



■id/3) 



12 Mi? 



+ 



— (-i? 2 + 12M 2 ) 
12 v ; 

M # %1 (V3) 



1 



M 



■ J-s/3 



1 %1 



(1/3) 



6 Mi? 



+ 



(-i? 2 + 12M 2 ^ 



V 



Mi?%l (1/3) 



%1 := i? 3 -54/3 M 2 i? 3 - 18 i? M 2 

+ 6V3V-M 2 H 2 + 16M 4 -H 6 p + 27H 6 (5 2 M 2 + 18H 4 (3 M 2 M 
Result through Jacobi sn — function 



u(0) = Ul + (u 2 - ui) JacobiSN(i ^2 a/m (u 3 - ui) + <$, J— — —f 



it 3 - ni 



As «i = ^- and U2 = — are the small values for the planets ( r\ and T2 are 
the perihelion and aphelion points) and u\ + U2 + u 3 = > we have 2M 
■U3 =1 and 

> u(theta) - u[l] =\ 

> (u[2] - u[l])*JacobiSN(l/2*theta+delta,0)~2; 



u(0) - m = (u 2 - ui) sin(- 9 + 8) 2 



u 2 +u 1 



U\ > 



that is the equation of orbital motion (see above) with e 
corresponds to the elliptical motion, u\ < corresponds to hyperbolical 
motion. The period of the orbital motion in the general case: 



20 



> (u[2]-u[l])/(u[3]-u[l]): 

> kernel := 2/sqrt ( (l-t~2) * (l-t~2*°/„) ) : #2 in the numerator\ 

> corresponds to sn~2-period 

> int (kernel, t=0..1); 

2EllipticK(./- U2 ~ Ul ) 



and can be found approximately for small u\ and u 2 as result of expansion: 

> series (2*EllipticK(x) , x=0,4): 

> convert (% , polynom) : 

> subs(x = sqrt( 2*M*(u[2]-u[l] ) ), I); 



7T+ -Tt M (u 2 - U\) 



From the expression for u{9 ) and obtained expression for the period of sn 2 
we have the change of angular coordinate over period: 

> %/(l/2*sqrt( 2*M*(u[3]-u[l] ) )): 

> simplify (%) ; 

1 7r(2 + Mu 2 -M«i)v / 2 

2 y/-M (-U3 + U1) 



But V2M(u 3 - m) = y/l -2M(ui-u 2 - ui) « 1 + M (2 u x + u 2 ) ■ And 
the result is 

> 2*Pi*(l + M*(u[2]-u[l])/2)*(l + M*(2*u[l]+u[2])) : 

> series('/„ u[l]=0,2): 

> convert polynom) : 

> series (%, u[2]=0,2) : 

> convert (%, polynom) : 

> expand (°/ ) : 

> expand (%-op (4//,)) : 

> f actor (°/ ) ; 



ir(2 + 3Mui+3Mu 2 ) 



The deviation of the period from 2 it causes the shift of the perihelion over 
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one rotation of the planet around massive star. 



> simplify (%-2*Pi) : 

> subs({u[l]=l/r[l] , u[2]=l/r[2]}, '/.) : 

> subs({r[l]=a*(l+e) , r [2]=a*(l-e)},G*%/c~2) :# we\ 

> returned the constants G and c 

> simplify (°/o) ; 

Git M 
a(l + e)(-l + e)c 2 

This result coincides with the expression, which was obtained on the basis 
of approximate solution of basic-equation. From the Kepler's low we can 
express this result through orbital period: 

> subs(M=4*Pi~2*a~3/T~2, '/„):# T is the orbital period 

> simplify ('/„); 

G^a 2 
_24 

T 2 (l + e)(-l + e) c 2 

2.6 Conclusion 

So, we found the expression for the Schwarzschild metric from the equiva- 
lence principle without introducing of the Einstein's equations. On this ba- 
sis and from Euler-Lagrange equations we obtained the main experimental 
consequences of GR-theory: the light ray deflection and planet's perihelion 
motion. 



3 Relativistic stars and black holes 
3.1 Introduction 

The most wonderful prediction of GR-theory is the existence of black holes, 
which are the objects with extremely strong gravitational field. The in- 
vestigation of these objects is the test of our understanding of space-time 
structure. We will base our consideration on the analytical approach that 
demands to consider only symmetrical space-times. But this restriction does 
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not decrease the significance of the obtained data because of the rich struc- 
ture of analytical results and possibilities of clear interpretation clarify the 
physical basis of the phenomenon in the strongly curved space-time. The 
basic principles can be found in [0, ||. 

3.2 Geometric units 

The very useful normalization in GR utilizes the so-called geometric units. 
Since the left-hand side of the Einstein equations describes the curvature 
tensor (its dimension is cm^ -2 ) ), the right-hand side is to have same dimen- 
sion. Let's the gravitational constant is G = 6.673 lO^ -8 ) = 1 and the 

° g s z 

light velocity is c = 2.998 10 10 = 1, 
then 

G/ c 2 = .74 25 10(- 28 ) cm/g = 1 

c 5 / G = 3.63 10 59 erg/s = 1 (power unit) 

G/c = 2.23 lO^ 18 -* Hz* cm 2 /g = 1 (characteristic of absorption) 
c 2 / ^G = 3.48 10 24 CGSE units (field strength) 
h/2 vr = 1.054 lot" 27 ) g* cm 2 /s = 2.612 10^ 66 ) cm 2 
elementary charge e = 1.381 10* -34 ) cm 
1 ps = 3.0856 10 18 cm 
1 eV = 1.324 10^ 61 ) cm 

There are the following extremal values of length, time, mass and density, 
which are useful in the context of the consederation of GR validity: 

= 1.616 10^ 33 ) cm (Planck length) 

= 5.391 10(" 44 ) s (Planck time) 
yjjjfe = 2.177 10(~ 5 ) g (Planck mass) 
2 -^4 = 5.157 10 93 gj cm 3 (Planck density) 

3.3 Relativistic star 

As stated above the first realistic metric was introduced by Schwarzschild 
for description of the spherically symmetric and static curved space. Let us 
introduce the spherically symmetric metric in the following form: 
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> restart: 

> with(tensor) : 

> with(plots) : 

> with(linalg) : 

> with(dif forms) : 

> coord := [t, r, theta, phi] :# spherical coordinates , \ 

> which will be designated in text as [0,1,2,3] 

> g_compts :=\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# metric components 

> g_compts [1 , 1] := -exp(2*Phi (r) ) : # component\ 

> of interval attached to d(t)~2 

> g_compts [2 , 2] := exp(2*Lambda(r) ) :# component\ 

> of interval attached to d(r)~2 

> g_compts [3,3] := r~2:# component of interval\ 

> attached to d(theta)~2 

> g_compts [4,4] := r~2*sin(theta) "1 : # component \ 

> of interval attached to d(phi)~2 

> g := create ( [-1,-1] , eval(g_compts) ) ;# covariant\ 

> metric tensor 

> ginv := invert( g, 'detg' );# contravariant\ 

> metric tensor 



table( [compts 



table( [compts 



- e (2*(r)) 







index-char = [—1, —1]]) 
ginv := 



e (2*(r)) 












e (2A(r)) o 








1 



e (2A(r)) 












r 2 

2 „;„/a\2 







o r 2 sm(ey _ 








1 



r 2 sin(6>) 2 



index-char = [1, 1]]) 
Now we can use the standard Maple procedure for Einstein tensor definition 
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> # intermediate values 

> Dig := dlmetric( g, coord ): 

> D2g := d2metric( Dig, coord ): 

> Cfl := Christoffell ( Dig ) : 

> RMN := Riemann( ginv, D2g, Cf 1 ) : 

> RICCI := Ricci( ginv, RMN ): 

> RS := Ricciscalar( ginv, RICCI ): 

> Estn := Einstein( g, RICCI, RS ):# Einstein tensor 

> displayGR (Einstein, Estn) ;# Its nonzero components 



The Einstein Tensor 



non — zero components : 



Gil = - 



e (2*(r)) ( 2 (9_ A ( r )) r + e (2 A(r)) _ 

r 2 e( 2A W) 



G22 



2(|-$(r))r-e( 2A ( r )) + 1 



G33 



;(2AW) 



l Hr) ~ 1 A(r) + r ^ *(r) + r (£ $(r)) 2 - r f $(r) & A(r) 



G44 = 



sin(e) 2 r 
' e (2A M ) 



I *(r) - & A(r) + r ^ *(r) + r (£ $(r)) 2 - r £ $(r) £ A(r) 



character : [—1,-1] 



In the beginning we will consider the star in the form of drop of liquid. In 
this case the energy- momentum tensor is T Mj „ = ( p + p ) u u + p g^ v 
(all components of u except for uq are equal to zero, and -1 = </ ' °) no uq ', 
the signature (-2) results in u a u a =1 and T Mi v = (p + p) u u - p g^v)'- 
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> T_compts :=\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# energy-momentum\ 

> tensor for drop of liquid 

> T_compts[l,l] := exp(2*Phi(r))*rho(r) : 

> T_compts [2,2] : = exp(2*Lambda(r) ) *p(r) : 

> T_compts [3,3] := p(r)*r~2: 

> T_compts [4,4] := p(r) *r~2*sin(theta) ~2 : 

> T := create ( [-1 , -1] , eval(T_compts)) ; 



T : = 



table( [compts 



- e (2-fM) p ( r ) 








p(r) r 





.2 



p(r) r 2 sin(#) 2 









index-char 



[-1 



, "I]]) 



To write the Einstein equations (sign corresponds to [||]) 



let's extract the tensor components: 



> Energy _momentum := get_compts (T) ; 

> Einstein := get_compts (Estn) ; 



Energy ^momentum := 



e (2*W) p ( r ) 

e( 2A W)p(r) 









p(r)r 2 







p(r) r 2 sin(^) 2 
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Einstein : = 

- e (2$(r)) ( 2 (|.A(r))r + %!-!) 
r 2 %l 

: o,- 2( *» (r)) r* 1+1 ,o,o 



0, 0, 



,0,0,0 



r ((* *W) " (* AW) + r *(r)) + r (£ $(r)) 2 - r (& <&(r)) (& A(r))) 



%1 



, 



0, 0, 0, 



sin(fl) 2 r «jL $(r)) - (& A(r)) + r S(r)) + r (& <l>(r)) 2 - r g(r)) (£ A(r))) 

%1 



%1 :=e (2A(r)) 

First Einstein equation for (0,0) - component is: 

> 8*Pi*Energy_momentum [1 , 1] + Einstein [1 , 1] : 

> expand (%/exp (Phi (r)) ~2) : 

> eql := simplify (%) = 0;#first Einstein equation 



eql := — - 



-8vrp(r) r 2 + 2e(~ 2A ( r )) {§- A(r)) r + 1 - e (~ 2A ( r )) 



= 



This equation can be rewritten as: 

> eql :=\ 

> -8*Pi*rho(r)*r~2+Diff (r* (l-exp(-2*Lambda(r) ) ) ,r) 



= 0; 



eql :=-8^/9(r)r 2 + (^-r(l-e ( - 2A(r)) )) =0 



The formal substitution results in: 
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> eql_n := subs(\ 

> r*(l-exp(-2*Lambda(r) ) ) = 2*m(r),\ 

> lhs(eql)) = 0; 



eql _n := -8 vr p(r) r 2 + ( — (2 m(r))) = 



> dsolve(eql_n,m(r)) ; 

m(r) = Jatt p(r)r 2 dr + _C1 



So, m(r) is the mass inside sphere with radius r. To estimate _Ci we have 
to express A from m: 

> r*(l-exp(-2*Lambda(r) ) ) = 2*m(r) : 

> expand (solve (%, exp(-2*Lambda(r) ) )); 

2m(r) 



Hence (see expression for g^ tU through A ) the Lorenzian metric in the 
absence of matter is possible only if _C1=0. 
Second Einstein equation for (l,l)-component is: 

> eq2 := simplify ( 8*Pi*Energy_momentum[2,2] +\ 

> Einstein [2 , 2] ) = 0;#second Einstein equation 

8tt e( 2A ( r » p(r) r 2 - 2 $(r)) r + e^ 2 ^ - 1 _ 

> eq2_2 := numer(\ 

> lhs(\ 

> simplify(\ 

> subs(exp(2*Lambda(r)) = l/(l-2*m(r)/r) ,eq2)))) = 0; 

eq2-2 := -8vr r 3 p(r) + 2 (— $(r)) r 2 - 4 (— $(r)) r m(r) -2m(r) =0 



As result we have: 

> eq2_3 := Dif f (Phi (r) ,r) = solve(eq2_2, dif f (Phi (r) ,r) ); 
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d T . . 47rr 3 p(r) + m(r) 
eq2.3 := — $(r) = -- 



Or r (— r + 2 m(r)) 

We can see, that the gradient of the gravitational potential $ is greater than 
in the Newtonian case $ = ™ , that is the pressure in GR is the source 
of gravitation. 

For further analysis we have to define the relativist equation of hydrody- 
namics (relativist Euler's equation): 

> compts := array ( [u_t ,u_r ,u_th,u_ph] ) : 

> u := create([l], compts) :# 4-velocity 

> Cf2 := Christoffel2 ( ginv, Cfl ): 

> (rho(r)+p(r))*get_compts(\ 

> cov_diff( u, coord, Cf2 ) ) [1 ,2] / (u_t) =\ 

> -dif f (p(r) ,r) ;# radial component of Euler equation, \ 

> u_r=u_th=u_ph=0 

> eq3 := Diff (Phi(r) ,r) = solve (*/., dif f (Phi (r) ,r) ); 

(P(r)+p(r))(| : *(r)) = -(| : p(r)) 



eg3 *(>■)-- ^ P(0 



dr p(r) + p(r) 

As result we obtain the so-called Oppenheimer-Volkoff equation for hydro- 
static equilibrium of star: 

> Diff(p(r),r) = factor( solve (rhs (eq3) =\ 

> rhs(eq2_3),diff(p(r),r))); 

d . (47rr 3 p(r) + m(r)) (p(r) + p(r)) 
dr ^ r (— r + 2m(r)) 



One can see that the pressure gradient is greater than in the classical limit ( 
^P = - ) an d this gradient is increased by pressure growth (numerator) 
and r decrease (denominator) due to approach to star center. So, one can 
conclude that in our model the gravitation is stronger than in Newtonian 
case. 
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Out of star m(r)=M, p=0 (M is the full mass of star). Then 

> diff (Phi(r) ,r) = subs ({m(r)=M,p(r)=0} ,rhs (eq2_3) ) ; 

> eq4 := dsolve(°/ , Phi(r)); 



— <3?(r) = 

dr v ' r(-r + 2M) 



eq4 :=$(r) = ln(r) + ^ ln(r - 2 M) + _CJ 



The boundary condition 



> = limit (rhs (eq4) ,r=infinity) ; 



= .CI 



results in 



> subs(_Cl=0,eq4) ; 



$(r) = -- ln(r) + ^ ln(r - 2 M) 



So, Schwarzschild metric out of star is: 



g_matrix := get_compts (g) ; 

_ e (2<&(r)) 

qjmatrix := ^ 








e (2A(r)) o 














r 2 sin(6>) 5 
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> coord := [t, r, theta, phi]: 

> sch_compts :=\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# metric components 

> sch_compts [1 , 1] := expand(\ 

> subs(Phi(r)=-l/2*ln(r)+l/2*ln(r-2*M) A 

> g_matrix [1 , 1] ) ):# coefficient of d(t)~2 in interval 

> sch_compts [2 , 2] := expand(\ 

> subs(Lambda(r)=-ln(l-2*M/r)/2,\ 

> g_matrix [2 , 2] ) ):# coefficient of d(r)~2 in interval 

> sch_compts [3,3] := g_matrix [3 , 3] :# coefficient\ 

> of d(theta)~2 in interval 

> sch_compts [4,4] := g_matrix[4,4] :# coefficient\ 

> of d(phi)~2 in interval 

> sen :=\ 

> create ( [-1 , -1] , eval (sch_compts) ) ; # Schwarzschild metric 



sch 



table( [compts 



-1 + 








2M 




1 



2M 

r 
















r z 
r 2 sin(#) 2 



index_char = [—1, —1]]) 



Now we consider the star, which is composed of an incompressible substance 
p = p = const (later we will use also the following approximation: p=( 7 
-1) p , where 7 =1 (dust), 4/3 (noncoherent radiation), 2 (hard matter) ). 
Then the mass is 

> m(r) = int(4*Pi*rhoO*r~2,r) ; 

> M = subs(r=R,rhs(%)) ;# full mass 



m(r) = - 7r pO r 



M = - vr pO R 3 
3 H 



and the pressure is 
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> diff (p(r) ,r) = subs(\ 

>{rho (r) =rhoO ,m (r) =4/3*Pi*rhoO*r~3} , \ 

> - (4*Pi*r ~3*p (r) +m(r) ) * (rho (r) +p (r) ) A 

> (r* (r-2*m(r) ) ) );# Oppenheimer-Volkof f equation 

> eq5 := dsolve(°/ 0) p(r)) ; 



(4 7r r 3 p(r) + - 7r pO r') (pO + p(r)) 
o 

. 8 rj. 

r yr — - it pOr 3 ) 



eq5 :=p(r)=pO( _ 3 + 8?rp0r 2 _ 9%1 ^ 

— 3 + 8 7r pO r z — 9 %1 

%1 . = e (-16_CJ7r P 0) 

> solve (\ 

> simplify(subs(r=R,rhs(eq5 [1] )))=0,\ 

> exp(-16*_Cl*Pi*rhoO) ) ; #boundary condition 

> solve (\ 

> simplify(subs(r=R,rhs(eq5 [2] )))=0,\ 

> exp(-16*_Cl*Pi*rhoO) ) ;#boundary condition 



-3 + 87rp0 R 2 



-3 + 87rp0 R 2 

> soil := simplify(subs(exp(-16*_Cl*Pi*rhoO) = 

> -3+8*Pi*rhoO*R~2 , rhs (eq5 [1] ) ) ) ; 

> sol2 := simplify(subs(exp(-16*_Cl*Pi*rhoO) = 

> -3+8*Pi*rhoO*R~2 , rhs (eq5 [2] ) ) ) ; 

soil := — 

4 



-3 + 12^p0i? 2 + y/9 - 24 7T pO R 2 - 24 vr pO r 2 + 64 tt 2 pO 2 r 2 i? 2 - 4 tt pO r 2 

-3-7rp0r 2 + 97TP0.R 2 

:= 4- 
4 

3- UnpOR 2 + ^9 - 24 tt pO i? 2 - 24 vr pO r 2 + 64 vr 2 pO 2 r 2 i? 2 + 4 vr pO r 2 

-3-7rp0r 2 + 97TP0.R 2 

In the center of star we have 
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> simplify (subs (r=0 , soli) ) ; 

> simplify (subs (r=0 , sol2) ) ; 

1 P (-3 + 12 7T pO R 2 + y/9 - 24 tt /?0 i? 2 ) 
~12 -l + STrpOi? 2 

1 p0(— 3 + 12TrpOR 2 -y/9- 24irpOR 2 ) 
~12 -l + 37rp0i? 2 

The pressure is infinity when the radius is equal to 

> R_critl =\ 

> simplify(solve(expand(denom(7o)/12)=0,R) [1] A 

> radical , symbolic) ; 

> R_crit2 =\ 

> simplify(solve(expand(denom(7 7 )/12)=0,R) [2] A 

> radical , symbolic) ; 



R-critl 



R_crit2 

> plot3d({subs (R=l , soil) , subs (R=l , sol2) } , \ 

> r=0. . I,rho0=0. .0. l,axes=boxed,title= < pressure < ) ;#only\ 

> positive solution has a physical sense 



1 V3^/F^0 
3 

1 y/3 

"3 ^pO 
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To imagine the space-time geometry it is necessary to consider an equator 
section of the star ( 9 = ir /2) at fixed time moment in 3-dimensional flat 
space. The corresponding procedure is named as "embedding". From the 

metric tensor, the 2-dimensional line element is ds 2 = - %„ ' - 2 



and for Euclidean 3-space we have ds 
investigate the 2-surface z=z(r). As dz 
line element: ds 2 



= [1+ (tr) 2 



dz' 

dz 
dr 



+ r dphi 

We will 
dr, one can obtain Euclidian 



dr 2 + r 2 dphi' 



dr 2 + r 2 dphi 2 



> subs (theta=Pi/2,get_compts (sen)) ; #Schwarzschild metric\ 

> dr-2*%[2,2] + dphi "1*% [3 ,3] = \ 

> (1+diff (z(r) ,r)"2)*dr"2 + r~2*dphi~2;#equality of\ 

> intervals of flat and embedded spaces 

> diff(z(r),r) = solve('/.,diff(z(r),r))[l]; 

> dsolve(%,z(r)) ;# embedding 



■1 + 








2M 




1 



1 - 



2M 



















r 2 sin(-7r) 2 
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dr 2 



+ r 2 dphi 2 



(1 + (— z(r)) 2 ) dr 2 + r 2 dphi 2 



2M 



1 - 



r 




V2 y/(r - 2M)M 
r-2M 



z(r) = 2 ^2rM-4M 2 + _C1 



Now let's take the penultimate equation and to express M through po = 
const: 

> rho_sol :=\ 

> solve(M = 4/3*Pi*rhoO*R"3,rhoO) ;#density from full mass 

> funl :=\ 

> Int(l/sqrt(r/((4/3)*Pi*rhoO*r~3)-l) ,r) ;# inside space 

> fun2 :=\ 

> Int(l/sqrt(r/((4/3)*Pi*rhoO)-l) ,r) ;# outside space 



rho_sol := - 



3 M 



4 ttR 3 





> fun3 := value (subs (rhoO=rho_sol, funl)) ;# inside 

> fun4 := value (subs (rho0=rho_sol,fun2)) ;# outside 



funS := 



B? + r 2 M 
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fun4 := 




-AM 



The resulting embedding for equatorial and vertical sections is presented be- 
low (Newtonian case corresponds to horizontal surface, that is the asymptote 
for r-> oo ). The outside space lies out of outside ring representing star 
border. 

> figl :=\ 

> plot3d(\ 

> subs(\ 

> M=l, subs (R=2.66*M, subs (r=sqrt(x~2+y~2) ,fun3))) ,\ 

> x=-5. . 5,y=-5. . 5,\ 

> grid= [100 , 100] , style=PATCHCONTDUR) :# the inside\ 

> of star [r in units of M] 

> fig2 :=\ 

> plot3d(\ 

> subs (M=l, subs (R=2.66*M, subs (r=sqrt (x~2+y~2) ,fun4))) ,\ 

> x=-5 . . 5 , y=-5 . . 5 , style=PATCHCONTDUR) : # the outside\ 

> of star 

> display (f igl , f ig2 , axes=boxed) ; 



We can see, that the change of radial coordinate dr versus the variation of 
interval dl= , dr 2M (dl is the length of segment of curve on the depicted 



1 



-1 

-2 
-3 
-4 
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surface) in vicinity of the star is smaller in compare with Newtonian case 
(plane z=0, where dr=dl). The star "rolls" the space so that an observer 
moves away from the star, but the radial coordinate increases slowly. But 
what is a hole in the vicinity of the center of surface, which illustrates the 
outside space? What happens if the star radius becomes equal or less than 
the radius of this hole R=2M? Such unique object is the so-called black hole 
(see below). 

But before consideration of the main features of black hole let us consider 
the motion of probe particle in space with Schwarzschild metric. We will 
base our consideration on the relativistic low of motion: g( a ' ^ p a pp + fi 2 
= (p is 4-momentum, fi is the rest mass). Then (if p =[ —E, r, 0, L ], 
where L corresponds to angle momentum, A is the affine parameter, second 
term describes the radial velocity): 

> -E"2/(l-2*M/r(tau))\ 

> + diff (r(tau) ,tau)~2/(l-2*M/r(tau)) +\ 

> L~2/r(tau) "2+1 = 0;# tau=lambda*mu, E=E/mu, L=L/mu 

> eq6 := dif f (r (tau) ,tau) ~2 =\ 

> expand(solve(°/ ,diff (r(tau) ,tau) ~2)) ;#integral of motion 

> V := sqrt(-factor(op(2,rhs(eq6)) + op(3,rhs(eq6)) +\ 

> op(4,rhs(eq6) ) + op(5,rhs(eq6)))) ;# effective potential 

E 2 (#r(r)) 2 L 2 

u K9r y " -| 1-1 = 

2M + i 2M + r (r) 2 

r(r) r(r) 



eq6 := — r r 2 = E 2 + — - - — — + — — - 1 

or r(r) v(t) 2 r(r) d 




V 



The effective potential determining the orbital motion is (we vary the angu- 
lar momentum L and r): 

> subs ({M=l,r (tau) =x, L=y}, V) : 

> plot ({subs (y=3 , '/.) , subs (y=4 , °/„) , subs (y=6 , °/„) , 1} , \ 

> x=0. . 30 , axes=boxed, color= [red, green, blue , brown] ,\ 

> title= t potentials < ,\ 

> view=0 . 6 . . 1 . 4) ;#level E=l corresponds\ 

> to a particle, which was initially motionless 
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potentials 




6 10 1^ 20 25 30 



We considered the motion in the relativistic potential in the first part, but 
only in the framework of linear approximation (weak field) . As result of such 
motion the perihelion rotation and light ray deflection appear. But now, in 
the more general case, one can see a very interesting feature of the relativistic 
orbital motion: unlike Newtonian case there is the maximum V max with the 
following potential decrease as result of radial coordinate decrease. If V max 
>E > 1 the motion is infinite (it is analogue of the hyperbolical motion in 
Newtonian case). V max >E=1 corresponds to parabolic motion. When E 
lies in the potential hole or E = V, we have the finite motion. The motion 
with energy, which is equal to extremal values of potential, corresponds to 
circle orbits (the stable motion for potential minimum, unstable one for the 
maximum). The existence of the extremes is defined by expressions: 

> numer( simplify ( diff( subs(r(tau)=r,V) , r) ) ) =\ 

> 0:# zeros of potential's derivative 

> solve (°/ ,r) ; 

1 (L + VL 2 - 12 M 2 ) L 1 (L - y/L 2 - 12 M 2 ) L 

2 M ' 2 M 



As consequence, the circle motion is possible if L > 2M , when r > 
3(2M). The minimal distance, which allows the circular unstable motion, is 
defined by 

> limit (1/2* (L-sqrt(L~2-12*M-2))*L/M, L=inf inity) ; 
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3M 



Lack of extremes (red curve) or transition to r < 3M results in the falling 
motion (the similar case is E > V max ). This is the so-called gravitational 
capture and has no analogy in Newtonian mechanics of point-like masses. 
The finite orbital motion in the case of the right local minimum existence 
(blue curve) is the analog of one in the Newtonian case, but has no elliptical 
character (see first part). 

As the particular case, let's consider the radial (£=0) fall of particle. The 
proper time of falling particle is defined by next expression 

> tau = Int(\ 1/sqrt (subs ({L=0,r (tau)=r} ,rhs (eq6) ) ) ,r) ; 

> tau = Int(l/(sqrt(2*M/r-2*M/R)) ,r) ;#R=2*M/(1-E~2) -\ 

> apoastr, i.e. radius of zero velocity 

> limit ( value( rhs(%) ) ,r=2*M ): 

> simplify (°/„) ; #So , this is convergent integral\ 

> for r— >2*M 




2M R +R R (- 2 ^V-(-^M + R)M + Rln(2) 
-R]n(-R + AM + 2 v / 2 v / -(- 2M + R) M)) y/2 



I ^-{-2M + R)M 
For the remote exterior observer we have -§f r = §t r '§f^ = '§i r JZ1K = 

r 

E r G (here ro is the time coordinate relatively infinitely remote observer) : 

> r[o] := Int(l/(l-2*M/r) ,r) = int(l/(l-2*M/r) ,r) ; 

> limit ( value ( rhs(%) ) ,r=2*M );#this is\ 

> divergent integral for r — >2*M 
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dr = r + 2 M ln(r — 2 M) 



r 



— signum(M) oo 



These equations will be utilized below. 



3.4 Degeneracy stars and gravitational collapse 

The previously obtained expressions for the time of radial fall have an in- 
teresting peculiarity if the radius of star is less or equal to 2M ( 2 ^2 M in 
dimensional case, that is the so-called "gravitational radius" R g ). The fi- 
nite proper time r of fall corresponds to the infinite "remote" time r Q , when 
r->2M. This means, that for the remote observer the fall does not finish 
never. It results from the relativistic time's slowing-down. The consequence 
of this phenomena is the infinite red shift of "falling" photons because of the 
red shift in Schwarzschild metric is defined by the next frequencies relation: 

I A_2M 

^ = 4^ = \ 90 ' o( f 2 \ = V 2 (here A r are the proper time intervals 

^2 An y 3o,o(n) L 2 m v ^ ^ 

between light flash in different radial points). Simultaneously, the escape 

velocity is equal to \J 2G R M that is the velocity of light for sphere with R= 
R g . So, we cannot receive any information from interior of black hole. 

How can appear the object with the radius, which is less than R g ? In the 
beginning we consider the pressure free radial (i.e. angular part of metric is 
equal to 0) collapse of dust sphere with mass M. 
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> r := 'r' : 

> E := >E> : 

> subs( r=r(t) ,get_compts(sch) ): 

> d(s)~2 =\ 

> °/„[l,l]*d(t)"2 + %[2,2]*d(r)~2;#Schwarzschild metric 

> -d(tau)~2 = collect(\ 

> subs( d(r)=diff (r(t) ,t)*d(t) ,rhs(%)) ,d(t)) ;#tau is\ 

> the proper time for the observer on the surface\ 

> of sphere 

> %/d(tau)~2; 

> subs({d(t)=E/(l-2*M/r(t)) ,d(tau)=l}, /,) ;#we used \ 

> d(t)/d(tau) = E/(l-2*M/r) 

> pot_l :=\ 

> factor( solve (°/„, (diff(r(t) ,t))~2) ) ; #"potential" 



2 M „ , 2 d(r) 2 



r(t) 



2M 



-d(r) 5 



-1 + 



2M (f r(t)) 2 
r(t) + 2M 
r(t)/ 



d(t) 2 



-1 + 



r(t) + 



1 - 



2M 



d(t) s 



d(r) 



-1 + 



-1 = 



2M (f r(t)) 2 
r(t) + 2M 
r(i)/ 



E 2 



(I-—) 2 
1 r(t) j 



_ (r(t) - 2 M) 2 (-r(t) + -E 2 r(t) + 2 M) 
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> plot({subs({E=0.5,M=l,r(t)=r},pot_l) ,0*r},\ 

> r=1.7. .3,axes=boxed,title= < (dr/dt)"2 vs r'); 

(dr/dt) A 2 vs r 



0.04 : 
0.02^ 

-0.02 : 




1.8 2 2.2 ^.4 2.6 2.8 



The collapse is the motion from the right point of zero velocity ^ r (this 
is the velocity for remote observer) to the left point of zero velocity. These 
points are 

> solve(subs(r(t)=r,pot_l)=0,r) ; 

2M, 2 .M. -2 ^— 

-1 + E 2 

that are the apoastr (see above) and the gravitational radius. Hence (from 
the value for apoastr and pot_l ) 

„ m 2 M \2 1 2 M U n2) /i 2M\2 ( 2M 2 M \ D / -, R g\ 

\dt') — E 2 ~ i — lM. — V 1 r > 1 1 i-Ha.)' 

The slowing-down of the collapse for remote observer in vicinity of i? 9 results 
from the time slowing-down (see above). But the collapsing observer has 
the proper velocity r : 

> pot_2 := simplify(pot_l*(E/(l-2*M/r(t)))~2) ;# 

> d/d(t)=(d/d(tau))*(l-2*M/r)/E 
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pot. 2 := 



-r(t) +E 2 x{t) +2M 



> plot({subs({E=0.5,M=l,r(t)=r},pot_2) ,0*r},\ 

> r=l .7. .3,axes=boxed,title =< (dr/dtau) ~2 vs r c ); 



(dr/dtau) A 2 vs r 




As result, the collapsing surface crosses the gravitational radius at finite 
time moment and lim r ^# 9 r =1 (that is the velocity of light), when E = 
1 (R= oo ). Time of collapse for falling observer is 

> Int(l/sqrt(R/r-l) ,r)/sqrt(l-E"2) ;#from pot_2,\ 

> R is apoastr simplify( value(%), radical, symbolic); 



1 



R 



dr 



„ rn — 7 „ .1 R-2r . 

i2 y7(i^7) + ^arctan(- ^_ r) ) 

2 y/l^W 



It is obvious, that the obtained integral is ^j=p J r R -7^ 1 = J[ R E i 



43 



TC M 



(1-E 2 ) ( 2> 



It is natural, that the real stars are not clouds of dust and the equilibrium is 
supported by nuclear reactions. But when the nuclear reactions in the star 
finish, the "fall" of the star's surface can be prevented only by pressure of 
electron or baryons. The equilibrium state corresponds to the minimum of 
the net-energy composed of gravitational energy and thermal kinetic 

energy of composition. Let consider the hydrogen ball. When the temper- 
ature tends to zero, the kinetic energy of electrons is not equal to zero due 
to quantum degeneracy. In this case one electron occupies the cell with vol- 
ume ~ A 3 , where A = 2 ^ p is the Compton wavelength ( p e is the electron's 
momentum). For non-relativistic electron: 



> E[e] := p [e] "2/m [e] ; #kinetic energy of electron 

> E[k] := simplify ( n [e] *R~3*subs ( p[e]=\ 

> h*n[e] ~(l/3)/(2*pi) ,E[e] ) );#full kinetic energy, \ 

> n[e] is the number density of electrons 

> E[k] := subs( n [e] =M/R~3/m [p] ,%) ;# m[p] is the mass\ 

> of proton. We supposed m[e]<<m[p] and n[p]=n[e] 

> E[g] := -G*M~2/R:# gravitational energy 

> E : = E [g] + E [k] ; # full energy 



E e :~- 



1 n e ( 5 /3) B 3 h 2 

Ek := 7 2 

4 7r z m e 



1 K B?m p 

tin 



)(5/3) R Z h 2 



J k ■= T 



7r^ nip 



I(J^_)(5/3 )jR 3 ft 2 

GM 2 4 { R 3 m v J 

E := p + ~2 

H 7r z m e 

> pot := -1/R+1/R~2:# the dependence of energy on R 

> plot(pot,R=0.5. . 10, \ 

> title= ' energy of degenerated star'); 
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1.5 
1 

0.5 




energy of degenerated star 



4 R 6 



8 



10 



One can see that the dependence of energy on radius has the minimum and, 
as consequence, there is the equilibrium state. 

> solve ( diff(E, R) = 0, R);# minimum of energy\ 

> corresponding to equilibrium state 



%lh 2 



1 



'I 



llV3%l\ 



2 Mm,, 



+ 



Mm r 



2 Mm p 2 ir 2 m e G' 2 



7r 2 m p m e G 



1 %1 



'^iVs%i\ 



2 Mm,, 



+ 



Mm r 



V 



7r 2 m p m e G 



%1 := (M 2 m p )W 
So, we have the equilibrium state with R m i n 



—5— and the star, 

A?) 



which exists in such equilibrium state as result of the termination of nuclear 
reactions with the following radius decrease, is the white dwarf . The value 
of number density in this state is 

> simplify ( subs(\ 

> R[min] =h~2/(G*M~ (1/3) *m[e] *m [p] "(5/3)) ,\ 

> n [e] =M/R [min] ~3/m [p] )); #equilibrium number density 
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M 2 G 3 m e 3 nip 4 



So, the mass's increase decreases the equilibrium radius (unlike model of 
liquid drop, see above) and to increase the number density (the maximal 
density of solid or liquid defined by atomic packing is ~20 g/ cm 3 , for white 
dwarf this value is ~ 10 7 g/ cm 3 !). But, in compliance with principle of 
uncertainty, the last causes the electron momentum increase ( p e n e ^~3^ 
h). Our nonrelativistic approximation implies p e « m e c or 

> simplify( rhs (%) " (1/3) *h ) - c*m[e] ;#this must be\ 

> large negative value 

> solve (%=0,M) ; 

,M 2 G 3 m e 3 m p \ nffl , 

( -) (1/3) h - cm e 



VGchhc VGchhc 
G 2 m p 2 G 2 m p 2 



The last result gives the nonrelativistic criterion M « ' c > . Therefore 

m p 2 G { ?> 

in the massive white dwarf the electrons have to be the relativistic particles: 

> E[e] := h*c*n [e] ~ (1/3) ; #E=p*c for relativistic particle 

> E[k] := simplify( n[e] *R~3*E [e] ); 

> E[k] := subs( n[e] =M/R~3/m[p] ,'/.) ; 

> E[g] := -G*M~2/R: 

> E : = E [g] + E[k];#this is dependence -a/R+b/R 

> #equilibrium state: 

> simplify( diff(E, R) , radical ): 

> numer (%) ; #there is not dependence on R therefore\ 

> there is not stable configuration with energy minimum 

> solve ( °/„=0, M ) ; 

E e := hcn e {1/3) 



E k ■= n e (4/3) R 3 hc 
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M 



E k := (J^-)^R 3 hc 



R y R^m p 



-(-GMm„ + (-=— ) (1/3) hcR) M 
it d m p 



0. 



VGchhc VGchhc 



G 2 m p 2 



G 2 m p 2 



We obtained the estimation for so-called critical mass M r 



- (fec) ( i ) 



=1.4 



-^sun (so-called Chandrasekhar limit for white dwarfs) . The smaller masses 
produce the white dwarf with non-relativistic electrons but larger masses 
causes collapse of star, which can not be prevented by pressure of degeneracy 
electrons. Such collapse can form a neutron star for 1.4 M sun < M < 3 M sun . 
Collapse for larger masses has to result in the black hole formation. 



3.5 Schwarzschild black hole 



Now we return to Schwarzschild metric. 



> get_compts (sen) ; 

2M 
-1 + 












1 



2M 

r 
















r 2 sin(fl) 2 



One can see two singularities: r=2M and r=0. What is a sense of first 
singularity? When we cross the horizon, sh\^\ and s/12,2 (i-e. 50,0 an d 
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gi, i ) change signs. The space and time exchange the roles! The fall gets 
inevitable as the time flowing. As consequence, when particle or signal cross 
the gravitational radius, they cannot escape the falling on r=0. This fact 

can be illustrated by infinite value of acceleration on r=2M, which is - - — — 

J ' so, o 

( T are the Christoffel symbols): 

> Dlsch := dlmetric( sen, coord ): 

> Cfl := Christoffell ( Dlsch ): 

> displayGR(Christoffell,°/„) ; 

The Christoffel Symbols of the First Kind 
non — zero components : 

M 

[22,2] = ~(r-2My 
[23, S]=r 
[24,4}=rsm(9) 2 
[33, 2] = -r 
[34,4] = r 2 sin(0) cos(6>) 
[44,2] =-rsm(9) 2 
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[44 , 3} = -r 2 sin(0) cos(#) 

> -get_compts(Cf 1) [1 , 1 ,2] /get_compts(sch) [1,1] ;#radial\ 

> component of acceleration 

M 



r l (-1 H / 



Such particles and signals will be expelled from the cause-effect chain of 
universe. Therefore an imaginary surface r= R g is named " event horizon ". 

The absence of true physical singularity for r= R g can be illustrated in the 
following way. The invariant of Riemann tensor R a , 13,^,6 i?( a '^' 7 ' 5 ) is 

> schinv := invert( sen, 'detg' ): 

> D2sch := d2metric( Dlsch, coord ): 

> Cfl := Christoffell ( Dlsch ): 

> RMN := Riemann( schinv, D2sch, Cfl ): 

> raise (schinv, RMN, 1) :#raise of indexes\ 

> in Riemann tensor 

> raise (schinv, %, 2) : 

> raise (schinv, °/ , 3) : 

> RMNinv := raise (schinv, %, 4) : 

> prod (RMN, RMNinv, [1,1] , [2,2] , [3,3] , [4,4]) ; 

M 2 

t&ble([compts = 48 — g-, index -char = []]) 



and has no singularity on horizon. This is only coordinate singularity and 
its sense is the lack of rigid coordinates inside horizon. True physical sin- 
gularity is r=0 and has the character of the so-called space-like singularity 
(the inavitable singularity for the observer crossing horizon). 

Now try to embed the instant equator section of curved space into flat space 
(see above): 

> z(r) [1] = int(sqrt(2*r*M-4*M~2)/(-r+2*M) ,r) ; 

> z(r) [2] = int(-sqrt(2*r*M-4*M~2)/(-r+2*M) ,r) ; 



z ( r ) l = -2 v / 2rM-4M 2 
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z(r) 2 = 2 ^2rM-4M 2 



> plot3d ({subs ({M=l , r=sqrt (x~2+y~2) } , rhs (%%) ) , \ 

> subs ({M=l , r=sqrt (x~2+y~2) } , rhs (%) ) } , \ 

> x=-10. . 10,y=-10. . 10,axes=boxed,style=PATCHC0NT0UR,\ 

> grid=[100,100] ,title= ' Schwarzschild black hole') ; 




This is the black hole, which looks as a neck of battle (Einstein-Rosen 
bridge): the inside way to horizon and the outside way to horizon are the 
ways between different but identically asymptotically flat universes ("worm- 
hole" through 2-dimensional sphere with minimal radius 2M). But since the 
static geometry is not valid upon horizon (note, that here t->t+dt is not 
time translation), this scheme is not stable. 

Let's try to exclude the above-mentioned coordinate singularity on the hori- 
zon by means of coordinate change. For example, we consider the null 
geodesies ( ds 2 =0) in Schwarzschild space-time corresponding to radial mo- 
tion of photons: dt 2 = — dr n g ^ 2 • Hence (see above considered expression for 

ro ): 

> t = combine( int(-l/(l-R[g] /r) ,r) ); 

t = -r-Rg ln(r - R g ) 
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Then if v is the constant, which defines the radial coordinate for fixed t, we 
have 



> t = -r - R[g]*ln( abs(-r/R[g]+l) ) + v;#module\ 

> allows to extend the expression for r<R[g] 



RgH 



Rg 



+ v 



Differentiation of this equation with subsequent substitution of dt 2 in ex- 
pression for interval in Schwarzschild metric produce 

> def f orm (f =0 , wl=l , w2=l , w3=l , v=l , R [g] =0 , r=0 , v=0) ; 

> d(t)~2 = expand(\ 

> subs( d(R[g])=0,\ 

> d( -r-R[g]*ln( r/R[g]-l )+v ) )~2);# dif f erentials\ 

> in new coordinates 

> subs( d(t) ~2=rhs(°/„) , sch_compts [1 , 1] *d(t) ~2 ) +\ 

> sch_compts [2,2] *d(r) ~2 + sch_compts [3,3] *d(theta) ~2 +\ 

> sch_compts [4,4] *d(phi) ~2: 

> collect ( simplify ( subs(R[g]=2*M,%) ),\ 
>{d(r) ~2,d(v) ~2}) ;#new metric 

d{tf = d{v? _2_div)im + r2 ^ 2 



Rg {r-Rgf 



(r-2M)d(v) 2 ~r 3 d(6) 2 - r 3 d(0) 2 + r 3 d(0) 2 cos(g) 2 

-\- 2 Q.[T} 0.[Vj 



So, we have a new linear element (in the so-called Eddington - Finkelstein 
coordinates) ds 2 = - ( 1 - ^ ) d V 2 + 2dvdr + r 2 d(9) 2 ( A(9) is the spher- 
ical part). The corresponding metric has the regular character in all region 
of r (except for r=0). It should be noted, that the regularization was made 
by transition to "light coordinate" v therefore such coordinates can not be 
realized physically, but formally we continued analytically the coordinates 
to all r>0. We can see, that for future directed (i.e. dv 2 >0) null ( ds 2 
=0) or time like ( ds 2 <0) worldlines dr<0 for r< R g that corresponds to 
above-mentioned conclusion about inevitable fall on singularity. 

In the conclusion, we consider the problem of the deviation from spher- 
ical symmetry for static black hole. Such deviation can be described by 
the characteristic of quadrupole momentum q. Erez and Rosen found the 
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corresponding static metric with axial symmetry: 



> coord := [t, lambda, mu, phi]: 

> er_compts :=\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# metric components 

> er_compts [1 , 1] :=\ 

> -exp(2*psi) :# coefficient of d(t)~2 in interval 

> er_compts [2,2] :=\ 

> M~2*exp(2*gamma-2*psi) * (lambda~2-mu~2) A 

> (lambda~2-l) : # coefficient of d(lambda)~2 in interval 

> er_compts [3,3] :=\ 

> M~2*exp(2 *gamma- 2*psi) * (lambda~2-mu~2) A 

> (l-mu~2) :# coefficient of d(mu)~2 in interval 

> er_compts [4,4] :=\ 

> M~2*exp(-2*psi)*(lambda~2-l)*(l-mu~2) : #coef f icient\ 

> of d(phi)~2 in interval 

> er := create ( [-1,-1] , eval(er_compts)) ;# axially\ 

> symmetric metric 

er : = table( [compts = 

-e^ 2 ^ 

M 2 e (2 7 -2V0 (A 2 - jj?) 

A2 _\ 

M 2 e (2 7 -2^) (A 2 -^ 2 ) 

\ ^ 

1 — fl z 

M 2 e(- 2 ^)(A 2 -l)(l 

, index -char = [—1, —1] 
]) 



where 



> fl := psi = l/2*( 

> (1+q* (3*lambda~2-l) * (3*mu~2-l) /4) *ln( (lambda- 1) A 

> (lambda+1) )+3/2*q*lambda*(3*mu~2-l) ); 

> f2 := gamma =\ 

> 1/2* (l+q+q'2) *ln( (lambda~2-l) / (lambda~2-mu~2) ) \ 

> -3/2*q* (l-mu~2) * (lambda*ln( (lambda-1) / (lambda+1) ) +2) \ 

> +9/4*q~2* (l-mu~2) * ( (lambda~2+mu~2-l-9*lambda~2*mu~2) *\ 

> (lambda~2-l)/16*ln( (lambda-1)/ (lambda+1) )~2+\ 

> (lambda~2+7*mu~2-5/3-9*mu~2*lambda~2) *lambda*\ 

> ln( (lambda-1) / (lambda+1) ) /4+l/4*lambda~2* (l-9*mu~2) +\ 

> (mu~2-l/3) ); 

> f3 := lambda = r/M-1; 

> f4 := mu = cos(theta); 

fl ■■= 

= ^ (1 + ^(3 A 2 - 1) (3^ 2 - l))ln(A^l) + ^A (3^ 2 - 1) 
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f2 ■■= i = l(i + q + g 2 )H^- 2 )- 

^(l-^ 2 )(Aln(^l) + 2) + ^ 2 (l-^)( 
1 (A 2 + ^ - 1 - 9 A V) (A 2 - 1) ln(^) 2 
+ i(A 2 + 7^ 2 -^-9AV)Aln(^)+ 
lA 2 (l-9^ 2 ) + ^ 2 -i) 

r 

fS ■= A = — - 1 
J M 



f4 ■= n = cos(6>) 



In the case g=0 we have 

> get_compts (er) : 

> map2 (subs , {psi=rhs (f 1) , gamma=rhs(f 2)} ,°/«) : 

> map2 (subs , q=0 , %) : 

> map2 (subs , {lambda=rhs (f 3) ,mu=rhs(f 4)} ,°/ ) : 

> map (simplify ,°/ ) ; 



r-2M 









rM 2 



r-2M 





-l + cos(0) 2 








-r 2 (-l + cos(6») 5 



That is the Schwarzschild metric with regard to f3 and f4- 

Now let's find the horizon in the general case of nonzero quadrupole moment. 
For static field ( ^g a ,p =0) the corresponding condition is 30,0 =0- Then 

> map2(subs,psi=rhs(f 1) ,er) : 

> er2 := map2 (subs ,gamma=rhs (f 2) ,°/ ) : 

> get_compts (er2) [1,1]; 
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e ((l+l/4 g (3A 2 -l) (3 /1 2 -l))ln(^i)+3/2 ? A(3 M 2 -l)) 



That is r=2M ( A =1). But the result for invariant of curvature is: 

> erinv := invert ( er2, 'detg' ): 

> Dler := dlmetric( er2, coord ): 

> D2er := d2metric( Dler, coord ): 

> Cfl := Christoffell ( Dler ): 

> RMN := Riemann( erinv, D2er, Cf 1 ) : 

> raise (erinv, RMN, 1) :#raise of indexes\ 

> in Riemann tensor 

> raise (erinv, °/ , 2) : 

> raise (erinv, %, 3) : 

> RMNinv := raise (erinv, °/ , 4) : 

> prod (RMN, RMNinv, [1,1] , [2,2] , [3,3] , [4,4]) : 

> get_compts (%) : 

> series (%,q=0, 2) :#expansion on q 

> convert (% , polynom) : 

> res := simplify (%) : 

In the spherical case the result corresponds to above obtained: 

> factor( subs(q=0,res) ):#spherical symmetry 

> subs (lambda=l, %) ;# this is 48*M~2/r~6 

3 1 

4 



There is no singularity. But for nonzero q (let's choose [i = for sake of 
simplification): 

> subs (mu=0, res) : 

> #L'Hospital's rule for calculation of limit 

> limit (diff (numer(%) , lambda) ,lambda=l) ; 

> limit (diff (denom(°/J , lambda) ,lambda=l) ; 

— signum((/) oo 







Hence, there is a true singularity on horizon that can be regarded as the 
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demonstration of the impossibility of static axially symmetric black hole 
with nonzero quadrupole momentum. Such momentum will be "taken away" 
by the gravitational waves in the process of the black hole formation. 



3.6 Reissner-Nordstrom black hole (charged black hole) 

The generalization of Schwarzschild metric on the case of spherically sym- 
metric vacuum solution of bounded Einstein-Maxwell equations results in 

> coord := [t, r, theta, phi]: 

> rn_compts := \ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# metric components 

> rn_compts [1 , 1] :=\ 

> -(l-2*M/r+Q~2/r~2) :# coefficient of d(t)~2 in interval 

> rn_compts [2 , 2] :=\ 

> l/(l-2*M/r+Q~2/r"2) :# coefficient of d(r)~2 in interval 

> rn_compts [3,3] := \ 

> g_matrix[3,3] :# coefficient of d(theta)~2 in interval 

> rn_compts [4,4] :=\ 

> g_matrix [4,4] :# coefficient of d(phi)~2 in interval 

> rn :=\ 

> create ( [-1 , -1] , eval (rn_compts) ) ; # Reissner-Nordstrom\ 

> (RN) metric 



rn : = 



table( [compts 



index_char = [— 1, 
]) 



-1 + 



-1] 



2M Q 2 










1 



2M Q 2 
1 + ^ 












r 2 

r 2 sin(6') 2 



Here Q is the electric charge. The metric has three singularities: r=0 and 

> denom(get_compts (rn) [2,2] ) = 0; 

> r_p := solve(%, r) [1] ; 

> r_n := solve (%%, r) [2] ; 



r 2 -2rM + Q 2 = 
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r_p :=M + ^M 2 -Q 2 



r_n := M — JM 2 - Q 



Let's calculate the invariant of curvature: 

> rninv := invert ( rn, 'detg' ): 

> Dlrn := dlmetric( rn, coord ): 

> D2rn := d2metric( Dlrn, coord ): 

> Cfl := Christoffell ( Dlrn ): 

> RMN := Riemann( rninv, D2rn, Cf 1 ) : 

> raise (rninv, RMN, 1) : #raise of indexes\ 

> in Riemann tensor 

> raise (rninv, %, 2) : 

> raise (rninv, %, 3) : 

> RMNinv := raise (rninv, %, 4) : 

> prod (RMN, RMNinv, [1,1] , [2,2] , [3,3] , [4,4] ) ; 

tabled^ = 8 S^-12,M<y + 7(y iniex _ ch ar = []]) 

> solve (get_compts (%) ,r) ; #nonphysical roots for\ 

> numerator with nonzero Q 

{l+ l -I^)Q 2 (l_I/V6)Q 2 

M ' M 



That is the situation like to one in Schwarzschild metric and two last sin- 
gularities have a coordinate character. Now let us plot the signs of two first 
terms in linear element (we plot inverse value for second term in order to 
escape the divergence due to coordinate singularities). 

> plot ({subs ({M=l,Q=l/2},get_compts(rn) [1,1]) ,\ 

> subs ({M=l , Q=l/2} , l/get_compts (rn) [2 , 2] ) } , \ 

> r=0.1..2, title =< signs of first and second terms of\ 

> linear element'); 
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signs of first and second terms of linear element 




-6-1 



One can see the radical difference from the Schwarzschild black hole. The 
space and time terms exchange the roles in region r_n < r < rjp (between 
so-called inner and outer horizons with go,o =0). But there are the usual 
signs in the vicinity of physical singularity, i. e. it has a time- like character 
and the falling observer can avoid this singularity. 

Next difference is the lack of coordinate singularities for M 2 < Q 2 . In this 
case we have the so-called naked singularity. One can demonstrate that there 
is no such singularity as result of usual collapse of charged shell with mass 
M and charge Q. The total energy (in Newtonian limit but with correction 
in framework of special relativity, M_0 is the rest mass) is 

> en := M(r) = M_0 + Q~2/r - M(r)~2/r;#we use\ 

> the geometric units of charge so that the Coulomb low is\ 

> G*Q_l*Q_2/r~2 

> sol := solve (en, M(r) ) ; 

Q 2 M(r) 2 

en := M(r) = M_0 + - — 

r r 

sol := --r + - yV 2 + 4M_0r + 4Q 2 , -- r - - ^ r 2 + 4 M _0 r + 4Q 2 
The choice of the solution is defined by the correct asymptotic lim r ^oo sol 
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= M_0: 



> if limit (sol [1] ,r=inf inity)=M_0 then true_sol :=\ 

> sol[l] fi: 

> if limit (sol [2] ,r=inf inity)=M_0 then true_sol :=\ 

> sol [2] fi: 

> true_sol; 

-\ r + \\/ r2 + AM - 0r + A Q 2 

> diff (true_sol,r) : 

> subs (M_0=solve (en , M_0) ,'/„): 

> simplify (%, radical, symbolic) ; 

-Q 2 + M(r) 2 
(r + 2M(r))r 

So, I M(r) 

~~ r(r+2M{r)) ■ c °U a P se i s possible if M decreases with 

decreasing R (domination of gravity over the Coulomb interaction in the 
process of collapse) that is possible, when M 2 > Q 2 . It should be noted, 
that the limit 

> limit (true_sol, r=0) ; 




resolves the problem of the infinite proper energy of charged particle. 

Now we will consider the pressure free collapse of charged sphere of dust by 
analogy with Schwarzschild metric. 



58 



> r := 'r' : 

> E := >E> : 

> subs( r=r(t) ,get_compts(rn) ): 

> d(s)~2 = 7„[l,l]*d(t)~2 + °/ [2,2]*d(r)"2;#RN metric 

> -d(tau)~2 =\ 

> collect(subs( d(r)=diff (r(t) ,t)*d(t) ,\ 

> rhs(°/ )) ,d(t)) ;#tau is the proper time for the observer\ 

> on the surface of sphere 

> %/d(tau)~2; 

> subs({d(t)=E/(l-2*M/r(t)+Q~2/r-2) A 

> d(tau)=l},°/„) ;#we used d(t)/d(tau)=E/(l-2*M/r+Q-2/r~2) 

> pot_l :=\ 

> factor( solve(°/„, (diff (r(t) ,t))~2) ) : #"potential"\ 

> for remote observer 

> pot_2 :=\ 

> simplify(pot_l*(E/(l-2*M/r(t)+Q'-2/r-2))"2) : #"potential"\ 

> for collapsing observer\ 

> d/d(t)=(d/d(tau))*(l-2*M/r+Q~2/r-2)/E 



d( S ) 2 = (-1 + 



2M Q 2 



r(t) T(t)- 



r)d(t) 2 + 



d(r) 



1 - 



2M 



Q 

r(t) ' Ht) 



+ 



-d(r) 2 



-1 + 



2M Q 2 



(f r(t)) 2 

2M Q 2 
~^t) + 7(t) 2 ) 



-1 + 



2M Q 2 

7(tj ~ r(ty< 



+ 



(|r(t)) 2 



2M Q 



+ 



d(t) 2 



-1 = 



d(T)S 



r(t) r(t)V 



-1 + 



2M 



Q2 



+ 



2M 
1 — + 



Q 



E 2 



-1 = 



r(t) r(t) 2 / 



1 r(t) ^ r 2 J 
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> plot({subs({E=0.5,M=l,Q=l/2,r(t)=r},pot_2) ,\ 

> 0*r},r=0. 11. .3,axes=boxed,title= < (dr/dtau)~2 vs r for\ 

> collapsing observer'); 

(dr/dtau) A 2 vs r for collapsing observer 




The difference from Schwarzschild collapse is obvious: the observer crosses 
the outer and inner horizons but does not reach the singularity because of 
the collapsar explodes as white hole due to repulsion with consequent rec- 
ollapse and so on. 

And at last, we consider the "extreme" case M 2 = Q 2 . 

> subs(Q~2=M~2, get_compts(rn)) ; 

> factor (% [1,1]); 



-1 + 



2M M 2 



r 











1 



2M M 2 
1 + 


















r 2 
r 2 sin(#) 2 



(r - Mf 
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So, we have one coordinate singularity in r=M. What happen with second 
horizon? Let's find the distance between horizons for fixed t and angular 
coordinates for RN-metric: 



> get_compts (rn) : 

> d(s)~2 = %[2,2]*d(r)~2;#RN metric 

> # or 

> d(s)~2 = d(r) ~2/expand( (l-r_p/r) * (l-r_n/r) );#second\ 

> representation of expression 



d( S ) 



d(rf 



2M Q 2 

1 I 7T 



d(r) 



1 2M Q< 

-L "To 



Hence, when r_p->r_n ( M 2 — Q 2 ->0) 

> r_p := 'r_p' : 

> r_n := 'r_n' : 

> s = Int(l/sqrt((l-r_p/r)*(l-r_n/r)) ,r=r_n. .r_p) ; 

> simplify ( value (rhs(°/ )) .radical, symbolic ); 




— r_nln(r_p — r_n) + — r_p\n(r_p — r_n) — 

- r_n ln( — rjp + r_n) — - r_pln( — rjp + r_n) 

we have s-> oo . So, there is the infinitely long Einstein-Rosen bridge 
(charged string) between horizons that means a lack of wormhole between 
asymptotically flat universes. This fact can be illustrated by means of em- 
bedding of equatorial section of static RN space in flat Euclidian space (see 
above). 
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> d(r)~2/(l-l/r)~2 =\ 

> (1+diff (z(r) ,r)~2)*d(r)~2;#equality of radial elements\ 

> of intervals, M=l 

> diff(z(r),r) = solve (%, dif f (z(r) ,r)) [1] ; 

> dsolve(%,z(r)) ;# embedding 



^V = (l + (^^)) 2 )d(r) 2 
(1 - -) 2 dr 



d V2 
— z(r) = — 
or r 



z(r) = 2 V2r - 1 + ln(V2r - 1 - 1) - ln(\/2r - 1 + 1) + _Ci 

> plot3d(subs(r=sqrt(x~2+y"2) ,\ 

> 2*sqrt(2*r-l)+ln(sqrt(2*r-l)-l)-ln(sqrt(2*r-l)+l)) ,\ 

> x=-10. . 10,y=-10. . 10,axes=boxed,style=PATCHC0NTDUR,\ 

> grid= [100, 100] ,title=' extreme RN black hole');\ 

> #we expressed arctanh through In 




The asymptotic behavior of RN-metric as r-> oo is Minkowski. For inves- 
tigation of the situation r->M let introduce the new coordinate (see, for 
example, [[HI ): 
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> rn_assym := subs(\ 

>{r=M*(l+lambda) ,Q~2=M~2},get_compts(rn) ) ; 



rn_assym : = 




A + l (A + l) 2 



2 1 















1 







2 1 











1 - 



A + l (A + l) 2 










M 2 (A + 1) 2 






M 2 (A + 1) 2 



sin(0) 2 



> ml :=\ 

> series (rn_assym [1 , 1] , lambda=0 ,3) ; # we keep only\ 

> leading term in lambda 

> m2 :=\ 

> series (rn_assym [2 , 2] ,lambda=0,3) ;# we keep only\ 

> leading term in lambda 

> d(s)~2 = convert (ml ,polynom) *d(t) ~2\ 

> + M~2*convert (m2 ,polynom) *d(lambda) ~2 +\ 

> M~2*d(0mega) ~2 ; #d(0mega) is spherical part 



This is the Robinson-Bcrtotti metric . The last term describes two - dimen- 
sional sphere with radius M (these dimensions are compactified in the vicin- 
ity of horizon) and the first terms corresponds to anti-de Sitter space-time 
(see below) with constant negative curvature. 

3.7 Kerr black hole (rotating black hole) 

In the general form the stationary rotating black hole is described by so- 
called Kerr-Newman metric, which in the Boyer-Linquist coordinates can 
be presented as: 



ml := -A 2 + 0(A 3 ) 



m2 ■= A~ 2 + 0(A _1 ) 



d(s) 2 = -X 2 d{t) 2 + 



M 2 d(A) 2 
A 2 



+ M 2 d(0) 2 
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> coord := [t, r, theta, phi]: 

> kn_compts :=\ 

> array (sparse , 1 .. 4, 1 .. 4) :# metric components , \ 

> a=J/M, J is the angular momentum 

> kn_compts [1 , 1] :=\ 

> -(Delta-a~2*sin(theta) ~2) /Sigma: #coefficient of d(t)~2 

> kn_compts [1 ,4] :=\ 

> -2*a*sin(theta) "2* (r~2+a~2-Delta) /\ 

> Sigma: # coefficient of d(t)*dphi 

> kn_compts [2 , 2] :=\ 

> Sigma/Delta: # coefficient of d(r)~2 

> kn_compts [3,3] :=\ 

> Sigma: # coefficient of d(theta)~2 

> kn_compts [4,4] :=\ 

> ( ( (r~2+a~2) ~2-Delta*a~2*sin(theta) ~2) /\ 

> Sigma) *sin(theta) ~2 : #coef f icient of d(phi)~2 

> kn := create ( [-1,-1] , eval (kn_compts) ) ; # Kerr-Newman\ 

> (KN) metric 

> #where 

> sub_l := Sigma = r~2+a~2*cos (theta) "1 ; 

> sub_2 := Delta = r~2-2*M*r+a~2+sqrt (Q~2+P"2) ;# P is\ 

> the magnetic (monopole) charge 

asm(9) 2 (r 2 + a 2 - A) 
E 




((r 2 + a 2 ) 2 - A a 2 sin(g) 2 ) sin(g) 2 



sub.l :=S = r 2 + a 2 cos(6') 2 
sub_2 := A = r 2 -2rM + a 2 + ^ Q 2 + P 2 

In the absence of charges, this results in Kerr metric. The obvious singular- 
ities are (except for an usual singularity of spherical coordinates 6 =0): 



kn := table ([compts = 



A-a 2 sm{9) 2 
S 











S 
A 











E 











index -char = [— 1, 


-1] 
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> r_p :=\ 

> solve( subs({Q=0,P=0},\ 

> rhs(sub_2))=0,r ) [1] ;#outer horizon 

> r_n :=\ 

> solve ( subs({Q=0,P=0} ,\ 

> rhs(sub_2))=0,r ) [2] ;#inner horizon 

> solve( subs({Q=0,P=0},rhs(sub_l))=0,theta ) ; 



r_p := M + \/M 2 - a 2 



r_n := M — \/M 2 - a 2 



- TT — /arcsinh(-), - ir + 7arcsinh(— ) 



The last produces r=0, 9 = § . 

As it was in the case of charged static black hole, there are three different 
situations: M 2 < a 2 , M 2 = a 2 , M 2 > a 2 . 

Let us consider the signs of <7o,o , <7i,i , 94, a ( M 2 > 

> plot3d(subs ({M=l , a=l/2 ,P=0 , Q=0} , subs (\ 
>{Sigma=rhs (sub_l) ,Delta=rhs (sub_2) } , \ 

> get_compts (kn) [1,1])) ,r=0 . 1 . . 4,theta=0 . .Pi , color=red) : 

> plot3d (subs ({M=l , a=l/2 , P=0 , Q=0 , theta=2*Pi/3} , \ 

> subs({Sigma=rhs(sub_l) ,Delta=rhs (sub_2) } , \ 

> l/get_compts(kn) [2,2] ) ) ,\ 

> r=0 . 1 . . 4 , theta=0 . . Pi , color=green) : 

> display Vh, 11 ,\ 

> title =< signs of first and second diagonal elements\ 

> of metric ' , axes=boxed) ; 

> plot3d (subs ({M=l ,a=l/2,P=0,Q=0} , subs (\ 
>{Sigma=rhs(sub_l) ,Delta=rhs(sub_2)} ,\ 

> get_compts(kn) [4,4])) ,r=-0.01. .0.01,\ 

> theta=Pi/2-0 . 01 . . Pi/2+0 . 01 , color=blue , \ 

> title='four diagonal element of metric in vicinity of\ 

> singularity' ,axes=boxed) ; 
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1-58 0.01 



One can see, that the approach to r=0 in the line, which differs from 9 = ^, 
corresponds to usual signs of diagonal elements of metric, i.e. the observer 
crosses r=0 and comes into region r<0 without collision with singularity. 
But from the second picture we can see the change of the g% t 3 -sign for r<0. 
Now (ft is the time like coordinate. But this coordinate has circle character 
and, as consequence, we find oneself in the world with closed time lines. The 
approach to r=0 in the line of 9 = ^ produces the change of go, o -sign, i.e. 
we find a true singularity in this direction. These facts demonstrate that 
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the singularity in Kerr black hole has a more complicated character than in 
above considered black holes. 



The more careful consideration gives the following results: 

1) M 2 < a 2 . There exist no horizon (r_p and r_n are complex), but the 
singularity in r=0, 9 = ^ keeps. To remove the coordinate singularity in 9 
=0 we introduce the Kerr-Schild coordinates with linear element 

> macro (ts= 't* ' ) : 

> ks_le := d(s)-2 = -d(ts)~2 + d(x)~2 + d(y)~2\ 

> + d(z)~2 + (2*M*r~3/(r~4+a~2*z~2))*\ 

> ( (r* (x*d(x) +y*d(y) ) -a* (x*d (y) -y*d (x) ) ) A 

> (r~2+a~2)+z*d(z)/r+d(ts) )~2; 

> x + I*y = (r + I*a)*sin(theta)*\ 

> exp(I*(Int(l,phi)+Int(a/Delta,r))) ; 

> z = r*cos (theta) ; 

> ts = Int(l,t) + Int((r~2+a~2/Delta) ,r) - r; 

ksJe := d(s) 2 = -d(t*) 2 + d(x) 2 + d(y) 2 + d(z) 2 



2 M r 3 C ^ d ^ + V d ^ ~ a ^ X d ^ - V d ^ + | d(t*)) 
_l r 2 + a 2 r 

r 4 + a 2 z 2 

x + Iy={r + Ia) sin(0) d 1 (f 1 d ^+I t *)) 



z = r cos(9) 



2 



t* = J\dt + Jr 2 + ^dr -r 



which is reduced to Minkowski metric by M->0. 

> int(subs( {P=0,Q=0},\ 

> subs (Delta=rhs (sub_2) , a/Delta) ),r): 

> x+I*y=(r+I*a)*sin(theta)*exp(I*(phi+ /,)) > 



x + Iy = (r + la) sin(0) e \ \ // 
When r=0, 9 = S the singularity is the ring x 2 + y 2 = a 2 , z = 0. 
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2) M 2 > a 2 . As before we have the ring singularity, but there are the 
horizons rjp and r_n. As an additional feature of Kerr metric we note here 
the existence of coordinate singularity: 

> get_compts (kn) [1 , 1] =0 ; 

> subs( Delta=rhs( sub_2 ),numer( lhs(°/J ) ): 

> solve( subs({Q=0,P=0},°/J = 0,r ); 

A - a 2 sin(g) 2 _ 
E 

M + ^M 2 -a 2 + a 2 sin(0) 2 , M - yj ' M 2 - a 2 + a 2 sin(#) 2 



The crossing of ellipsoid r_l = M+-\/M 2 — a 2 cos(#) 2 produces the change of 
go, -sign. As it was for Schwarzschild black hole this fact demonstrates the 
lack of static coordinates under this surface, which is called "ergosphere" and 
the region between rjp and rA is the ergoregion. The absence of singularity 
for r r<r suggests that the nonstatic behavior results from entrainment of 
observer by black hole rotation, but not from fall on singularity, as it takes 
place for observer under horizon in Schwarzschild black hole. 

3.8 Conclusion 

So, the elementary analysis by means of basic Maple 6 functions allows to 
obtain the main results of black hole physics including conditions of collapsar 
formation, the space-time structure of static spherically symmetric charged 
and uncharged black holes and stationary axially symmetric black hole. 



4 Cosmological models 
4.1 Introduction 

I suppose that the most impressive achievements of the GR-theory belong 
to the cosmology. The description of the global structure of the universe 
changes our opinion about reality and radically widens the horizon of knowl- 
edge. The considerable progress in the observations and measures allows to 
say that we are living in a gold age of cosmology. Here we will consider some 



68 



basic conceptions of relativistic cosmology by means of analytical capabili- 
ties of Maple 6. 

4.2 Robertson- Walker metric 

We will base on the conception of foliation of 4-dimensional manifold on 
3-dimcnsional space-like hyperplanes with a time-dependent geometry. Now 
let's restrict oneself to the isotropic and homogeneous evolution models, i.e. 
the models without dependence of curvature on the observer's location or 
orientation. We will suppose also, that such foliation with isotropic and ho- 
mogeneous geometry and energy-momentum distribution exists at each time 
moment. These space-like foliations are the so-called hyperplanes of homogeneity. 
The time vectors are defined by the proper time course in each space point. 

Let us consider 3-geometry at fixed time moment. As it is known, in 3- 
dimensional space the curvature tensor has the following form: 

Riklm 

Ril 9km - Rim 9kl + Rkm 9il " Rkl 9im + f ( 9im 9kl " 9il 9km )■ 

But if (as it takes a place in our case) there is not the dependence of curvature 
on space-direction in arbitrary point, our space has a constant curvature 
{Schur's theorem) that results in 

Riklm = f ( 9im 9kl - 9il 9km ), (1) 

where R is the constant Ricci scalar. 

Now we will consider the so-called Robertson- Walker metric : 

9ij = Ft s ■ ■ h~i Ti (^) 

(1 1 ' J 8 ) 2 
For this metric the spatial curvature tensor is: 
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> restart: 

> with(tensor) : 

> with(linalg) : 

> with(dif forms) : 

> coord := [x, y, z] :# coordinates 

> g_compts :=\ 

> array (symmetric , sparse , 1 .. 3, 1 .. 3) :# metric components 

> g_compts [1 , 1] := l/\ 

> (l+R*(x~2+y~2+z~2)/8) ~2:# component of interval\ 

> attached to d(x)~2 

> g_compts [2 , 2] :=\ 

> l/(l+R*(x~2+y~2+z~2)/8) ~2:# component of interval\ 

> attached to d(y)~2 

> g_compts [3,3] :=\ 

> l/(l+R*(x~2+y~2+z~2)/8)~2:# component of\ 

> interval attached to d(z)~2 

> g := create ( [-1 , -1] , eval(g_compts) ) ;# covariant\ 

> metric tensor 

> ginv := invert( g, 'detg' ):# contravariant\ 

> metric tensor 

> Dig := dlmetric( g, coord ):# calculation of \ 

> curvature tensor 

> D2g := d2metric( Dig, coord ): 

> Cfl := Christoffell ( Dig ) : 

> RMN := Riemann( ginv, D2g, Cf 1 ) : 

> displayGR(Riemann, RMN); 



g := table( [index -char = [—1, —1], compts 



%1 
%1 







%1 



]) 



%1 



(l + ^i?(x 2 + y 2 + z 2 )) 2 



The Riemann Tensor 



non — zero components : 



R1212 = 2048 



R 



(8 + Rx 2 + Ry 2 + Rz 2 )* 



R1313 = 2048 



R 



{S + Rx 2 + Ry 2 + Rz 2 Y 
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R2323 = 2048 



(8 + Rx 2 + Ry 2 + Rz 2 ) i 
character : [—1, —1, —1, — 1} 



Now we try to perform the direct calculations from Eq. (1): 

> A := get_compts (g) : 

> B := array (1. .3, 1..3, 1..3, 1..3): 

> for i from 1 to 3 do 

> for j from 1 to 3 do 

> for k from 1 to 3 do 

> for 1 from 1 to 3 do 

> B[i,j,k,l] :=\ 

> eval(R*(A[i,k]*A[j ,1] - A [i , 1] *A [ j ,k])/2) :# Eq. l.\ 

> These calculations can be performed also by means of\ 

> tensor [prod] and tensor [permute_indices] 

> if B[i,j,k,l]<>0 then \ 

> print( [i,j ,k,l] ,B[i,j ,k,l] );# non-zero components 

> else 

> f i: 

> od: od: od: od: 



[1, 2,1,2],- r 



2 (l + ^R(x 2 +y 2 + Z 2 )Y 



[1,2,2,1],-! _ R 



2 {l + ^R(x 2 + y 2 + z 2 )) A 



[1,3,1,3] ' 



2 (l + \R{x 2 +y 2 + z 2 )Y 



[1, 3, 3, 1], ' ^ 



2 {l + \R{x 2 + y 2 + z 2 )Y 
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[2, 1, 1, 2], -\ j 

2 ( 1+ jR ( x 2 +y2 + z 2))4 



[2,1,2,1] ,± 1 U 



(1 + -R(x 2 +y 2 + z 2 )f 



[2, 3, 2, 3], 



R 



2 (1 + 1^2+^2 + ^)4 



[2,3,3,2], -1 T * 



(l + ^O^ + ^ + z 2 )) 4 



[3,1,1,3],— T 



(l + ^R(x 2 + y 2 + z 2 )) 4 



[3, 1, 3, 1], \ j 

2 ( 1 + ^( x 2 +y2 + z 2))4 



[3, 2, 2, 3], -1 * 

2 ( 1 + i jR ( x 2 +y2 + z 2))4 



[3, 2, 3, 2], \ 

2 (1+ 0(^.2+^ + ^)4 



The direct calculation on the basis of Eq. (1) and metric (2) results in the 
tensor with symmetries [i,j,k,l] = [k,l,i,j] = [j,i,k,T\ = -[ i,j,l,k] and [i,j,k,l] + 
[i,l,j,k] + [i,k,l,j] = 0. These symmetries and values of nonzero components 
correspond to curvature tensor, which was calculated previously from (2). As 
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consequence, metric (2) satisfies to condition of constant curvature (Schur's 
theorem, Eq. (1)). We will base our further calculations on this metric. 



The next step is the choice of the appropriate coordinates. We consider two 
types of coordinates: spherical and pseudospherical. 

In the case of the spherical coordinates the transition of Eq. (2) to a new 
basis results in: 

> g_sp :=\ 

> simplify (subs (\ 

>{x=r*cos (phi) *sin(theta) ,y=r*sin(phi) *sin(theta) ,\ 

> z=r*cos(theta)},get_compts(g))) :# This is transf ormation\ 

> of coordinates in metric tensor 

> Sp_compts := array (1 .. 3 , 1 .. 3) :# matrix of\ 

> conversion from (NB!) spherical coordinates 

> Sp_compts [1 , 1] := cos (phi) *sin(theta) : 

> Sp_compts [2 , 1] := sin(phi) *sin(theta) : 

> Sp_compts [3, 1] := cos(theta) : 

> Sp_compts [1 , 2] := r*cos (phi) *cos (theta) : 

> Sp_compts [2 , 2] := r*sin(phi) *cos (theta) : 

> Sp_compts [3,2] := -r*sin(theta) : 

> Sp_compts [1 ,3] := -r*sin(phi) *sin(theta) : 

> Sp_compts [2 , 3] := r*cos (phi) *sin(theta) : 

> Sp_compts [3,3] := 0: 

> Sp := eval(Sp_compts) : 

> simplify(\ 

> multiply (transpose (Sp) , g_sp, Sp) ) ; #transition\ 

> to new coordinates in metric 

u Wi 

64— 



%1 

-64 



r 2 (-l + cos(#) 2 ) 



%1 

%1 := (8 + R r 2 cos(c/>) 2 sin(0) 2 + R r 2 sin(0) 2 sm(9) 2 + Rr 2 cos(9) 2 ) 2 
The denominator is 

> denom(%[l,l] ) : 

> expand (%) : 

> simplifyC/ ,trig) : 

> factor (*/,); 

(i?r 2 + 8) 2 
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After substitution r* = r ^Jl^- (we will omit asterix and suppose R > 0), 
the metric is converted as: 

> coord := [r, theta, phi] :# spherical coordinates 

> g_compts_sp :=\ 

> array (symmetric , sparse , 1 .. 3, 1 .. 3) :# metric components 

> g_compts_sp [1 , 1] :=\ 

> (2/R)/(l+r~2/4) ~2:# component of interval\ 

> attached to d(r)~2 

> g_compts_sp [2 , 2] :=\ 

> (2/R)*r~2/(l+r~2/4) ~2:# component of interval\ 

> attached to d(theta)~2 

> g_compts_sp[3,3] :=\ 

> (2/R)*r~2*sin(theta)~2/(l+r~2/4)~2:# component of\ 

> interval attached to d(phi)~2 

> g_sp := create ( [-1 , -1] , eval (g_compts_sp) ) ;# covariant\ 

> metric tensor 



gsp := table ( [index -char = [—1, —1], 
1 - 



compts 



R(l + 



1 



,2\2 



R(l + -r 2 ) 2 
4 



r 2 sin(0) 2 

i?(l + ^r 2 ) 2 
4 



]) 



The denominator's form suggests the transformation of radial coordinate: 
sin( x ) = r /{l + r<2 /4). Then for dr 2 -component we have 

> diff (sin(chi) ,chi)*d(chi) =\ 

> simplify(diff (r/(l+r"2/4) ,r))*d(r) ; 

> lhs(°/ )~2 = rhs(°/„)-2; 

> subs (cos (chi) ~2=l-r~2/ (l+r~2/4) ~2 , lhs ) \ 

> = rhs(%) : 

> factor ('/„); 



cos(x) d(x) = -4 ■ 



-4 + r 2 ) d(r) 
(4 + r 2 ) 2 



cos( X ) 2 d( X ) 2 = 16. 



-4 + r 2 ) 2 d(r) 
(4 + r 2 ) 4 
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( r _2) 2 ( r + 2) 2 d(x) 2 
(4 + r 2 ) 2 



= 16 



(r-2) 2 (r + 2) 2 d(r) 2 
(4 + r 2 ) 4 



that leads to a new metric (<i % = 



(l+£) 2 



> coord := [r, theta, phi] :# spherical coordinates 

> g_compts_sp :=\ 

> array (symmetric , sparse , 1 .. 3, 1 .. 3) :# metric components 

> g_compts_sp [1 , 1] := \ 

> 2/R:# component of interval attached to d(chi)~2 

> g_compts_sp [2 , 2] :=\ 

> (2/R) *sin(chi) "1 : # component of interval\ 

> attached to d(theta)~2 

> g_compts_sp[3,3] :=\ 

> (2/R) *sin(chi) ~2*sin(theta) ~2 : # component of\ 

> interval attached to d(phi)~2 

> g_sp := \ 

> create ( [-1 , -1] , eval (g_compts_sp) ) ; # covariant\ 

> metric tensor 



Here R/2 is the Gaussian curvature K (K>0) so that the linear element is: 

> l_sp := d(s)~2 = a~2*(d(chi)~2 +\ 

> sin(chi) ~2* (d(theta) ~2+sin(theta) "2*d(phi) ~2) ) ; 



l_ sp : = d(s) 2 = a 2 (d( X ) 2 + sin( X ) 2 (d(9) 2 + S m(9) 2 d(0) 2 )) 



Here a = J -g is the scaling factor. Note, that substitution r = a sin(x) 



produces the alternative form of metric: 

> subs(\ 

>{sin(chi)=r/a,d(chi)~2=d(r)~2/(a~2*(l-r-2/a-2))},\ 

> l_sp):#we use d(r)=a*cos(chi)*d(chi) 

> expand (°/ ) ; 



gsp := t&ble([index_char = [—1, —1], compts = 
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1_ ^2 



For negative curvature: 



> coord := [r, theta, phi] :# spherical coordinates 

> g_compts_hyp :=\ 

> array (symmetric , sparse , 1 .. 3, 1 .. 3) :# metric components 

> g_compts_hyp [1 , 1] :=\ 

> (-2/R)/(l-r~2/4)~2:# component of interval\ 

> attached to d(r)~2 

> g_compts_hyp[2,2] :=\ 

> (-2/R)*r~2/(l-r~2/4)~2:# component of interval\ 

> attached to d(theta)~2 

> g_compts_hyp [3,3] :=\ 

> (-2/R)*r~2*sin(theta)~2/(l-r~2/4)~2:#component of\ 

> interval attached to d(phi)~2 

> g-hyp :=\ 

> create ( [-1,-1] , eval (g_compts_hyp) ) ;# covariant metric\ 

> tensor 



gJiyp := table( [index-char 
1 



[—1, —1], compts 



R{l-\r 2 ) 2 
4 



]) 







R(l- -r 2 )2 
4 







-2 





r 2 sin(6>) 2 



R(l 



1 



,2^2 



where R is the module of Ricci scalar and new radial coordinate is / r 




Let's try to use the substitution sh( X ) = r /{l- f 2 /4)- Then for dr 2 
-component we have 

> diff (sinh(chi) ,chi)*d(chi) =\ 

> simplify(diff (r/(l-r~2/4) ,r))*d(r) ; 

> lhs(7 )-2 = rhs(%)~2; 

> subs(cosh(chi)"2=l+r"2/(l-r-2/4)-2,lhs(7,))\ 

> = rhsC/J : 

> f actor (*/,); 
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cosh(x) d(x) = 4 



(4 + r 2 )d(r) 
(-4 + r 2 ) 2 



cosh(x) 2 d( X ) 2 = 16 



(4 + r 2 ) 2 d(r) 2 
(-4 + r 2 ) 4 



(4 + r 2 ) 2 d( X ) 2 
(r - 2) 2 (r + 2) : 




(4 + r 2 ) 2 d(r) 2 
(r - 2) 4 (r + 2) 4 



And finally 



> coord := [r, theta, phi] :# spherical coordinates 

> g_compts_hyp :=\ 

> array (symmetric , sparse , 1 .. 3, 1 .. 3) :# metric components 

> g_compts_hyp [1 , 1] :=\ 

> -2/R:# component of interval\ 

> attached to d(chi)~2 

> g_compts_hyp [2,2] :=\ 

> (-2/R) *sinh(chi) ~2 : # component of interval\ 

> attached to d(theta)~2 

> g_compts_hyp [3,3] :=\ 

> (-2/R) *sinh(chi) ~2*sin(theta) ~2 : # component\ 

> of interval attached to d(phi)~2 

> g-hyp :=\ 

> create ( [-1,-1] , eval (g_compts_hyp) ) ; # covariant\ 

> metric tensor 



Here R/2 is the Gaussian curvature K (K<0). So, the linear element is 

> l_hyp := d(s)~2 = a"2*(d(chi)"2 +\ 

> sinh(chi) ~2* (d(theta) ~2+sin(theta) ~2*d(phi) ~2) ) ; 



gJiyp := 

table([index-char = [—1, —1], compts = 




sinh(y) 2 

-2 

R 

2 sinh(x) 2 sin(g) 2 



]) 



R 
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Lhyp := d(s) 2 = a 2 (d( X ) 2 + sinh( X ) 2 (d(6) 2 + sin(0) 2 d(0) 2 )) 



where a = y is the imaginary scale factor. The substitution r = a sinh( 
X ) results in: 

> subs(\ 

>{sinh(chi)=r/a,d(chi)'-2=d(r)-2/(a"2*(l+r'"2/a-2))},\ 

> l_hyp):#we use d(r)=a*cosh(chi) *d(chi) 

> expand (°/ ) ; 

d(s )2 = + r 2 d{()) 2 + r 2 sin( ^ )2 d(0) 2 

1 + - 



At last, for flat space: 

> l_flat := d(s)"2 = a~2* (d(chi) "2 +\ 

> chi~2*(d(theta)~2+sin(theta)~2*d(phi)"2)) ; 



Lflat := d(s) 2 = a 2 (d( X ) 2 + X 2 (d(6) 2 + sm(6) 2 d(0) 2 )) 
Here a is the arbitrary scaling factor. 



To imagine the considered geometries we can use the usual technique of 
embedding. In the general case the imbedding of 3-space is not possible 
if the dimension of enveloping flat space is only 4 (we have 6 functions 
gij , but the number of free parameters in 4-dimensional space is only 4). 
However, the constant nonzero curvature allows to embed our curved space 
in 4-dimensional flat space as result of next transformation of coordinates 
(for spherical geometry): 

> defform(f=0,wl=l,w2=l,w3=l,v=l,chi=0,theta=0,phi=0) ; 

> e_l := Dw~2 = a~2*d(cos (chi) ) ~2 ; 

> e_2 := Dx~2 =\ 

> (a*d(sin(chi) *sin(theta) *cos (phi) ) ) "2; 

> e_3 := Dy~2 =\ 

> (a*d(sin(chi) *sin(theta) *sin(phi) ) ) ~2; 

> e_4 := Dz~2 =\ 

> (a*d(sin(chi)*cos(theta)))~2; 
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e_l := Dw 2 = a 2 sin(x) 2 d(x) 2 



e_2 := Dx 2 = a 2 

(sin(0) cos(^) cos(x) d(%) + sin(x) cos(<?!>) cos(0) d{9) — 
sin(x) sin(0) sin(0) d((f>)) 2 

e_3 := By 2 = a 2 

(sin(0) sm(<j>) cos(x) d(%) + sin(x) sin(</>) cos(6 l ) d(#)+ 
sin(x) sin(0) cos(0) d(</>)) 2 

:= L>z 2 = a 2 (cos(0) cos(x) d(x) - sin(x) sin(6>) d((9)) 2 

> #so we have for Euclidian 4-metric: 

> simplify(\ 

> rhs(e_l) + rhs(e_2) + rhs(e_3) + rhs(e_4), trig); 

a 2 (d(</>) 2 cos(x) 2 cos(fl) 2 + d(6) 2 + d(4>) 2 - d(6) 2 cos{ X ) 2 - 
d(0) 2 co S (0) 2 -d(0) 2 cos( X ) 2 + d(x) 2 ) 

> collect (°/.,d(phi) "2) ; 

a 2 (1 - cos(fl) 2 - cos(x) 2 + cos(x) 2 cos(fl) 2 ) d(0) 2 + 
a 2 (d(0) 2 +d( X ) 2 -d(0) 2 cos(x) 2 ) 

This expression can be easily converted in Lsp, i. e. the imbedding is correct. 
Since 

> w~2 + x~2 + y~2 + z"2 =\ 

> simplify(\ 

> (a*cos(chi))~2 +(a*sin(chi)*sin(theta)*cos(phi))~2 +\ 

> (a*sin(chi)*sin(theta)*sin(phi))~2 +\ 

> (a*sin(chi) *cos (theta) ) ~2) ; 

w 2 + x 2 + y 2 + z 2 = a 2 



the considered hypersurface is the 3-sphere in flat 4-space. The scale factor 
is the radius of the closed spherical universe. 

In the case Lhyp there is not the imbedding in flat Euclidian 4-space, but 
there is the imbedding in flat hyperbolical flat space with interval 
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-d w 2 + d x 2 + d y 2 + d z 2 



As result we have 

> w~2 - x~2 - y~2 - z"2 =\ 

> simplify(\ 

> (a*cosh(chi) ) ~2 - (a*sinh(chi) *sin(theta) *cos (phi) ) ~2 -\ 

> (a*sinh(chi) *sin(theta) *sin(phi) ) "2 -\ 

> (a*sinh(chi)*cos(theta))~2) ; 

2 2 2 2 2 

w — x — y — z = a 



This is expression for 3-hyperboloid in fiat 4-space. 



4.3 Standard models 

Now we turn to some specific isotropic and homogeneous models. The as- 
sumption of time-dependent character of the geometry results in the time- 
dependence of scaling factor aft). Then 4-metric (spherical geometry) is 

> coord := [t, chi, theta, phi]: 

> g_compts := array (symmetric , sparse , 1 . .4, 1 .. 4) : 

> g_compts[l,l] := -l:#dt"2 

> g_compts[2,2] := a(t) ~2 : #chi~2 

> g_compts [3,3] :=a(t)~2*sin(chi)~2:#dtheta~2 

> g_compts [4,4] :=a(t)~2*sin(chi)~2*sin(theta)~2:#dphi~2 

> g := create ( [-1 , -1] , eval(g_compts) ) ; 

> ginv := invert ( g, 'detg' ): 

g := table( [index-char = [—1, —1], compts = 

"-10 

a(t) 2 

a(t) 2 sin(x) 2 " 

a(t) 2 sin(x) 2 sin(#) 2 _ 

The Einstein tensor is 

> Dig := dlmetric( g, coord ): 

> D2g := d2metric( Dig, coord ): 

> Cfl := Christoffell ( Dig ) : 

> RMN := Riemann( ginv, D2g, Cf 1 ) : 

> RICCI := Ricci( ginv, RMN ): 

> RS := Ricciscalar( ginv, RICCI ): 

> G : = Einstein( g, RICCI, RS ); 
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G := table( [index_char = [— 1, —1], 
compts = 

-2 sin( X ) 2 -3(| a(t)) 2 sin( X ) 2 - 1 + cos( X ) 2 q q q 
a(t) 2 sin(x) 2 

n -2 a(t) sin( X ) 2 (g, a(t)) - (j a(t)) 2 sin( X ) 2 - 1 + cos( X ) 2 

sm(x) 2 

0,0,2 a(t) sin(x) 2 a(t)) + sin( X ) 2 + ( * a(t)) 2 sin( X ) 2 , 
.0,0,0,2 a(t) sin( X ) 2 sin(fl) 2 (|* a(t)) + sin(£) 2 sin( X ) 2 + (| a(t)) 2 sin( X ) 2 sin(0) 2 

]) 

Let's the matter is defined by the energy-momentum tensor T^ v = (p+ p ) 
Up u u +p v : 

> T_compts :=\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# energy-momentum\ 

> tensor for drop of liquid 

> T_compts[l,l] :=rho(t): 

> T_compts[2,2] := a(t)~2*p(t) : 

> T_compts [3,3] : = a(t) ~2*sin(chi) ~2*p(t) : 

> T_compts [4,4] :=\ 

> a(t) ~2*sin(chi) ~2*sin(theta) ~2*p(t) : 

> T := create ( [-1,-1] , eval(T_compts) ) ; 

T := table([index-char = [—1, —1], compts = 

' p(t) 

a(t) 2 p(t) 

a(t) 2 sin(x) 2 p(t) 

a(t) 2 sin(x) 2 sin((9) 2 p0) 

]) 

Then the first Einstein equation Go,o - A go, o = - 8 ir T^o (we add the 
so-called cosmological constant A , which can be considered as the energy 
density of vacuum): 
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> get_compts (G) : 

> 51 [1,1]: 

> el := simplify (°/ ) : 

> get_compts (g) : 

> 51 [1,1]: 

> e2 := simplify ('/„) : 

> get_compts (T) : 

> 51 [1,1]: 

> e3 := simplify (50 : 

> E[l] := expand(el/(-3)) = -8*Pi*expand(e3/(-3) ) +\ 

> expand (e2*Lambda/ (-3) ) ;#first Einstein equation 



1 (|a(i)) 2 8 . . 1 



Second equation is: 

> get_compts (G) : 

> y.[2,2]: 

> el := simplify (°/ ) : 

> get_compts (g) : 

> y.[2,2]: 

> e2 := simplify (%) : 

> get_compts (T) : 

> y.[2,2]: 

> e3 := simplify (%) : 

> E[2] :=\ 

> expand(el/a(t)~2) = -8*Pi*expand(e3/a(t)~2) +\ 

> expand(e2*Lambda/a(t) ~2) ;#second Einstein equation 

|U(t) (^a(t)) 2 1 
E 2 := 2 M^l + LotWj + = -Svrp t + A 
a(i) a(t) 2 a(i) 2 



Two other equations are tautological. 

After elementary transformation we have for the second equation: 
> expand(simplify(E[2]-E[l])/2) ; 



&(t) ^ w 3 3 



Similarly, for the hyperbolical geometry we have: 
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> coord := [t, chi, theta, phi]: 

> g_compts := array (symmetric , sparse , 1 . .4, 1 .. 4) : 

> g_compts[l,l] := -l:#dt~2 

> g_compts [2,2] := a(t) ~2 : #chi~2 

> g_compts [3,3] :=a(t)~2*sinh(chi)~2:#dtheta~2 

> g_compts [4,4] :=\ 

> a(t) ~2*sinh(chi) ~2*sin(theta) ~2 : #dphi~2 

> g := create ( [-1,-1] , eval(g_compts) ) ; 

> ginv := invert ( g, 'detg' ): 

g := table( [index_char = [—1, —1], compts = 

"-10 " 

a(i) 2 

a(t) 2 sinh(x) 2 
a(t) 2 sinh(x) 2 sm(6») 2 

]) 

> Dig := dlmetric( g, coord ): 

> D2g := d2metric( Dig, coord ): 

> Cfl := Christoffell ( Dig ) : 

> RMN := Riemann( ginv, D2g, Cf 1 ) : 

> RICCI := Ricci( ginv, RMN ): 

> RS := Ricciscalar( ginv, RICCI ): 

> G : = Einstein ( g, RICCI, RS ); 

:= table ( [index -char = [—1, —1], compts = 

-2 sinh( X ) 2 + 3 (| a(t)) 2 sinh( X ) 2 + 1 - cosh( X ) 2 q q q 
a(i) 2 sinh(x) 2 

2 a(Q sinh( X ) 2 (g a(t)) + (| a(t)) 2 sinh( X ) 2 + 1 - cosh( X ) 2 

smh(x) 

0,0,2 a(t) sinh(x) 2 (-gr a(t)) - sinh( X ) 2 + (| a(i)) 2 sinh( X ) 2 , 
0,0,0,2 a(t) sinh(x) 2 sin(fl) 2 ( J^r a(t)) - sin(#) 2 sinh(x) 2 + (f a(t)) 2 sinh(x) 2 sin(fl) 2 

> T_compts :=\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# energy-momentum\ 

> tensor for drop of liquid 

> T_compts[l,l] :=rho(t): 

> T_compts [2,2] : = a(t) ~2*p(t) : 

> T_compts[3,3] := a(t) ~2*sinh(chi) ~2*p(t) : 

> T_compts [4,4] : = a(t) ~2*sinh(chi) ~2*sin(theta) ~2*p(t) : 

> T := create ( [-1,-1] , eval(T_compts) ) ; 
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T := table([index _char = [— 1, —1], compts = 

p(t) 

a(t) 2 p(t) 

a(t) 2 sinh(x) 2 p(t) 

a(t) 2 sinh(x) 2 sin((9) 2 p(t) 

]) 

> get_compts (G) : 

> 51 [1,1]: 

> el := simplify (%) : 

> get_compts (g) : 

> 51 [1,1]: 

> e2 := simplify (50 : 

> get_compts (T) : 

> 31 [1,1]: 

> e3 := simplify (%) : 

> E[l] :=\ 

> expand(el/(-3)) = -8*Pi*expand(e3/ (-3) ) +\ 

> expand(e2*Lambda/ (-3) ) ; #f irst Einstein equation 



1 (4a(()) 2 8 ,, 1, 

> get_compts (G) : 

> 7.12,2]: 

> el := simplify C/o) : 

> get_compts (g) : 

> 7. [2, 2]: 

> e2 := simplify (*/,) : 

> get_compts(T) : 

> °/.[2,2]: 

> e3 := simplify (%) : 

> E[2] :=\ 

> expand(el/a(t)"2) = -8*Pi*expand(e3/a(t) ~2) +\ 

> expand(e2*Lambda/a(t) "2) : #second Einstein equation 

> expand (simplify (E [2] -E [1] )/2) ; 



a(t) 



The second equation is identical with one in spherical geometry. 

In general, the first and second differential equations for scaling factor in 
our ideal universe: 
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> basic_l :=\ 

> (diff (a(t) ,t)/a(t))-2 = -K/a(t)~2 + Lambda/3 +\ 

> 8*Pi*rho (t) /3 ; #f irst basic equation for\ 

> homogeneous universe 

> basic_2 :=\ 

> diff(a(t), ( $ ( (t,2))/a(t) =\ 

> -4*Pi*(p(t)+rho(t)/3)+l/3*Lambda;# second\ 

> basic equation for homogeneous universe 



7§z a(t) 1 1 

basic _ 2 : = 9L_L± = -4tt (p(t + - p(t)) + - A 
&(t) 3 3 



Here K=+l for spherical, -1 for hyperbolical and for flat geometries. The 
left-hand side of the first equation is the square of Hubble constant H. If 
we denote recent value of this parameter as Hq , then Hq = 65 s k ^ pc (ps is 
the parallax second corresponding to distance about of 3.26 light-years or 
3* 10 13 km) Q. 

When K= A =0, the universe is defined by the so-called critical density p c 
= ^Mnl =7.9* 10 (_30) — . Let's denote the recent relative values fi* = ^ 
of the matter, radiation, cosmological constant and curvature as: 

n M = , n* = §§g* ,fl A = jA, ,n K = - . 

Here ao is the recent scaling factor. At this moment the contribution of the 
radiation to density is negligible « £Im ■ To rewrite the evolutionary 
equation in terms of fij we have to take into account the time dependence 
of these parameters: A does not depend on a{t) (in framework of standard 
model), the dependence of curvature ~ \ , density of matter ~ \ , but for 
radiation ~ . The last results from the change of density of quanta ~ -\ 
and their energy ~ - . The resulting equation can be written as: 



> basic_3 := H(a)"2 =\ 

H[0]~2* (Omega [R]* a [0]~4/a"4 + Omega [M] *a[0] "3/a"3\ 
+ Omega [Lambda] + Omega [K] *a[0] ~2/a~2) ; 



basiC-3 := H(a) z = H ( 1 ~ h "A H 5 — ) 



a 4 a 3 1 a 2 



As result of these definitions we have the cosmic sum rule: 



85 



1 = Um + + + 



Let us transit to the new variables: y= , r = #0 ( * — *o )> where to is 
the present time moment. Then 

> basic_4 :=\ 

> diff (y(tau) ,tau) = sqrt(l + (l/y(tau)-l)*Omega[M] +\ 

> (y(tau) ~2-l) *0mega [Lambda] ); #we neglect the radiation\ 

> contribution and use the cosmic sum rule.\ 

> d(tau)=H[0]*d(t) — > a[0] *H[0] *(d(y)/d(tau) )/a =\ 

> (d(a)/d(t))/a = H 



basics := y(r) = Jl + (-^ - 1) Q M + (y(r) 2 - 1) H A 



Now we will consider some typical standard models. For review see, for 
example, [p^] , or on-line survey [13|. 



4.3.1 Einstein static 



> rhs (basic_l) = 0: 

> rhs(basic_2) = 0: 

> allvalues(\ 

> solve({%, °/o°/ } , {a(t) , Lambda}) );# static universe 



{A = 12vrp(t) + 4vrp(t), a(i) = J- --}, 

V 47rp(t) + 4tt p(t) 



{A = 12vrp(t) + 4vrp(t), a(t) = - J- --} 

y 47rp(t) + Att p(t) 

At this moment for the matter p=0 and the contribution of the radiation is 
small, therefore: 



a 2 = -n— , A = 4ir p 

4 7T p ' r 



This solution has a very simple form and requires K = +1 (the closed 
spherical universe) and positive A (repulsive action of vacuum) , but does not 
satisfy the observations of red shift in the spectra of extragalactic sources. 
The last is the well-known fact demonstrating the expansion of universe (the 
increase of scaling factor). 
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4.3.2 Einstein-de Sitter 



= 1 and, as consequence = = 0. This is the flat universe 
described by equation 

> y(tau)~(l/2)*diff (y(tau) ,tau) = 1; 

> dsolve({%, y(0)=l}, y(tau)): 

> allvalues (°/ ) ; 

^)(§- T y(r)) = i 

y(r) = 

\(12r + 8)( 2 / 3 ) , y(r) = (- \ (12 r + 8)( 1 /3) + * / ^ (12 T + 8 )(V3) )2 ? 

y(r) = 

(-^ (12 r + 8)( 1 /3) _ 1 / ^3 (12 T + 8 )(V3) )2 

> plot (1/4* (12*tau+8)~ (2/3) , \ 

> tau= -2/3.. 2, axes=B0XED,\ 

> title= 'Einstein-de Sitter universe'); 



Einstein-de Sitter universe 




We have the infinite expansion from the initial singularity (Big Bang). 
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From the value of scaling factor y(r) = - — T —^- — it is possible to find the 
universe's age: 

> =\ 

> (12*H[0]*(0-t[0]) + 8)~(2/3)/4:# y(0)=0 

> solvea,t[0]); 

2 1 
3lo 

> evalf( (2/3) *3*10~ 19/65/60/60/24/365) ;#[yr] 

.9756859072 10 10 



The obtained universe's age does not satisfy the observations of the star's 
age in older globular clusters. 



4.3.3 de Sitter and anti-de Sitter 

The empty universes ( fi^ =0) with R > (de Sitter) and R < (anti-de 
Sitter). Turning to basicA we have for de Sitter universe: 

> assume (lambda, positive) : 

> subs ({rho (t) =0 , K=l , Lambda=lambda*3} , \ 

> basic_l) : #lambda=Lambda/3 

> expand(%*a(t)~2) ; 

> dsolve(%,a(t)) ; 

(|a(t)) 2 = -l+a(t) 2 r 



a(i) = j_ , a(i) 



%2 2 . %1 2 , 

(1 + _^)%1 + 

a (t) = i %1 , a (t) - 



2 



2 ^02^ ' W 2 %1 VF 
%2 := e( tv/X ") 

It is obvious, that the first two solutions don't satisfy the basic-2. Another 
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solutions give: 



where C is the positive or negative constant defining the time moment 
responding to the minimal scaling factor. 



> plot (subs (lambda=0 . 5 , \ 

> cosh(sqrt (lambda) *tau) /2/ sqrt (lambda) ) , \ 

> tau=-3. . 3,title =t de Sitter universe axes=boxed) ; 



de Sitter universe 




-3-2-10 1 2 

tau 



For anti-de Sitter model we have: 

> assume (lambda, negative) : 

> subs ({rho (t) =0 , K=-l , Lambda=lambda*3} , \ 

> basic_l) : #lambda=Lambda/3 

> expand(°/ *a(t)~2) ; 

> dsolve(%,a(t)) ; 



(^a(i)) 2 = l + a(t) 2 A~ 
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\/-A~ \/-A 
a(i) = — — — , a(t) 



a(t) = 
a(t) = 



A~ ' v ' A~ ' 

sin(t y/^A^ - _CJ \/ z A") 

sin(t y/^A" - .CI V^y) 

7^ 



De Sitter and anti-de Sitter universes play a crucial role in many modern 
issues in GR and cosmology, in particular, in inflationary scenarios of the 
beginning of universe evolution. 



4.3.4 Closed Friedmann-Lemaitre 

VLm > 1 A = and, as consequence VLk < (K=+l). This is the closed 
spherical universe. 

> y(tau)*(diff (y(tau) ,tau)~2+b) =\ 

> Omega [M] ;# b=0mega[M]-l > 



y(r)((^y(r)) 2 + b) = n M 



> solve ( °/ ,diff (y(tau) ,tau) ); 



v/y(r) (-y(r) b + Q M ) y/jij) (-y(r) b + Sl M ) 



y(r) 



y(r) 



> d(t) =sqrt(y/b/ (Omega [M]/b - y))*d(y); 



d(t) 



H—-y) 



d(y) 



Introducing a new variable we have : 

> sqrt (y/ (Omega [M] /b - y)) = tan (phi) ; 

> y = solve (°/,, y) ; 

> y = convert (rhs (°/ ) , sincos) ; 
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i 



tan(</>) 



tan((/>) 2 Qm 
b (1 + tan(<£) 2 ) 



sin(c/>) 2 fi M 



cos(» 2 &(l + 



sin(c/)) 2 
cos(0) 2 



But this is y = " M 7 W2 



Qm (1-cos(2 0)) 
2b 



And 



> defform(f=0,wl=0,w2=0,v=l,phi=0,y=0) ; 

> d(y) = d( (Omega [M]/b)*sin (phi) ~2 ); 



d(y) = 2 sin(<£) cos(</>) (d(0) ^) + sin(» 2 d(^) 



Then our equation results in 

> subs(\ 

> y=0mega[M]*sin(phi)~2/b,sqrt(y/(b*(0mega[M]/b-y))) ) 

> simplify (%) ; 



/ -l + cos(» 2 
V 6cos(0) 2 

> d(t) = simplify(\ 

> (Omega [M] /b) *tan(phi) *2*sin(phi) *cos (phi) ) \ 

> *d(phi)/sqrt(b) ; 



d(t) = -2 



^ M (-l + cos(0) 2 )d(0) 

6(3/2) 



Hence, the age of universe is: 
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> #initial conditions: y(0)=0 — > phi=0 

> int( Omega [M] *(l-cos(2*x))/ (Omega [M]-l)~ (3/2) ,\ 

> x=0..phi): 

> sol_l := t = simplify (°/ ) ; #age\ 

> of universe vs phi 

1 ^ M (-2 + sin(2 0)) 

2 (Q M - 1)0/2) 

This equation in the combination with y( <\> ) is the parametric representation 
of cycloid: 

> y = 0mega[M]*(l-cos(2*phi))/(2*(0mega[M]-l)) : 

> plot( [subs (Omega [M] =5 ,rhs(sol_l)) ,\ 

> subs (Omega [M] =5, rhs(°/ )) ,phi=0. .Pi] , axes=boxed,\ 

> title= ' closed Friedmann-Lemaitre universe'); 



closed Friedmann-Lemaitre universe 




In this model the expansion from the singularity changes into contraction 
and terminates in singularity (Big Crunch). 

The age of universe is: 
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> #initial conditions: y(tO)=l — > phi=? 

> 1 = sin(phi) ~2*0mega[M]/ (Omega [M] -1) : 

> sol_2 := solve (%, phi) ; 

> subs(phi=sol_2 [1] ,rhs(sol_l)) ;#age \ 

> of universe vs Omega [M] 

> plotO/,, Omega [M] =1.01. . 10, axes=BOXED, \ 

> title='age of universe [1/H[0] ]'); 



soL2 := arcsin( 



-), -arcsm( — ) 



M 



M 



V^M (pM — 1) 



-. Qm (—2 arcsin( — ) + sin(2 arcsin( 

1 



. \/$ Im Pm — 1) 



))) 



(n M - 1)0/2) 
age of universe [ 1/H[0] ] 




4 6 
Omega[M] 



So, the increase of £Im decreases the age of universe, which is less than 
in Einstein-de Sitter model. The last and the modern observations of mi- 
crowave background suggesting Q,k ~0 do not agree with this model. 



4.3.5 Open Friedmann-Lemaitre 



VLm < 1 A = and, as consequence Six > (K= -1). This is the open 
hyperbolical spherical universe. From the above introduced equation we 
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have 



> y(tau)*(diff (y(tau) ,tau)~2 - a) =\ 

> Omega [M] ;# a = -b = 1-Omega [M] >0 



y(r)((^-j(r)f-a) = n M 



> solve ( %,dif f (y(tau) ,tau) ); 

v / y( r ) (y(t) a + ^m) Vy(r) (y(r) a + Q M ) 

> dif f (y(tau) , tau) =sqrt ( (Omega [M] +\ 

> a*y(tau))/y(tau)) ; 



* y(T) = V W) 

> dsolve(°/ , y(tau)): 

> sol_l := simplify( subs (a=l-Omega [M] ,%) ,\ 

> radical, symbolic) ;#implicit solution for y(tau) 

soLl := --(-2r%10 M + ^Mln(-- 



-n M - 2y(r) + 2y(r) Q M - 2 y/-y(r) (-y(r) + y(r) M - »m) %! ■ 

v/l — £7m 



-2_Ci %1 fi M + 2 r %1 - 2 V-y(r) (-y(r) + y(r) M - «m) %1+ 
2_Ci %1) /(V1-»M (fijif - 1)) = 



%1 := Vl-^M 

> #definition of _C1 

> subs (tau=l , sol_l) : 

> subs(y(l)=l,%) : 

> solve (% , _C1) : 

> simplify( subs(_Cl=%,sol_l) ): 

> subs (y (tau) =y, °/„) : 

> sol_2 := solve (°/ , tau) ;#implicit solution\ 

> with defined _C1 
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soL2 := -\ (- 4 %1 M + n M ln( _ (^M ~ 2 - 2 %1) 2 , 
4 \ Om - 1 



-In - 



(JIm + 2y - 2y n M + 2 v / -y(-y + y^M-^M)%i) z 



M 



n 



M 



+ 4 A /-y(-y + yfiM-fiM)%i /(Vi - n M (0 M - l)) 



%i := VI -»m 

This solution can be plotted in the following form: 

> subs (Omega [M] =0.3, sol_2) : 

> plotC/,, y=0.01..2, axes=B0XED, colour=RED,\ 

> title= 'hyperbolical universe\ 

> (time vs scaling factor) ' ) ; 

hyperbolical universe (time vs scaling factor) 




0.2 0.4 0.6 0.8 



The extreme case Qm ->0 (almost empty world) 
> limit (sol_2, Omega [M] =0) ; 




results in the linear low of universe expansion and gives the estimation for 
maximal age of universe in this model: 
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> evalf (3*10"19/65/60/60/24/365) ;#[yr] 



.146352886110 



The dependence on VLm looks as: 

> - (subs (y=0 , sol_2) -subs (y=l , sol_2) ) : #age\ 

> of universe [ 1/H[0] ] 

> plot(%, Omega [M] =0. . l,axes=boxed,\ 

> title='age of universe [ 1/H[0]] f ); 




It is natural, the growth of matter density decreases the universe age. If esti- 
mations of Q.m give 0.3 (that results in the age « 12 gigayears for considered 
model), then this model does not satisfies the observations. Moreover, it is 
possible, that the nonzero curvature is not appropriate (see below). 

4.3.6 Expanding spherical and recollapsing hyperbolical universes 

This model demonstrates the possibility of the infinite expansion of universe 
with spherical geometry in the presence of nonzero cosmological constant. 
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> with(DEtools) : 

> DEplot (subs ({Omega [M] = 1 , Omega [Lambda] =2 . 59} , \ 

> basic_4) , [y(tau)] ,tau=-2.5. .1, [[y(0)=l]] , stepsize=0 . 01 , \ 

> axes=boxed , linecolor=red, arrows=none , \ 

> title =< loitering expansion'); 

loitering expansion 



y(tau)2 : 




-2.5 -2 -1.5 -1 -0.5 0.5 1 

tau 



We selected the parameters, which are close to ones for static universe (J5m 
=1) but the cosmological term slightly dominates tt\ > 2Qm ■ There is the 
delay of expansion, which looks like transition to collapse, but the domina- 
tion of 17a causes the further expansion. This model can explain the large 
age of universe for large Hq and to explain the misalignment in the "red 
shift-distance" distribution for distant quasars. 

The next example demonstrate the collapsing behavior of an open (hyper- 
bolical) universe: 

> pi :=\ 

> DEplot (subs ({Omega [M] =0.5, Omega [Lambda] =-1} , \ 

> basic_4) , [y(tau)] ,tau=-0.65. .0.7, [[y(0)=l]] ,\ 

> stepsize=0 . 01 , axes=boxed, linecolor=red, \ 

> arrows=none) :#expansion 

> p2 :=\ 

> DEplot (subs ({Omega [M] =0.5, Omega [Lambda] =-1} , \ 

> lhs(basic_4)=-rhs(basic_4)) , [y(tau)] ,\ 

> tau=0.71. .1.98, [[y (0.71) =1.36]] , stepsize=0 . 01 , \ 

> axes=boxed,linecolor=red,arrows=none) :#we have to\ 

> change the sign in right-hand side of basic_4 and\ 

> to change the initial condition 

> with(plots) : 

> display(pl,p2,title =< recollapsing open universe'); 
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recollapsing open universe 

1.4 

V 

0.8 ~_ 
y(tau) ; 
0.6- 

0.4^ 

0.2 

0z 

-0.5 0.5. 1 1.5 

tau 



4.3.7 Bouncing model 

Main feature of the above considered nonstatic models (except for de Sitter 
model) is the presence of singularity at the beginning of expansion and, it 
is possible, at the end of evolution. On the one hand, the singularity is not 
desirable, but, on the other hand, the standard model of Big Bang demands 
a high-density and hot beginning of expansion. The last corresponds to the 
conditions in the vicinity of the singularity. The choice of the parameters 
allows the behavior without singularity: 

> pi :=\ 

> DEplot (subs ({Omega [M] =0.1, Omega [Lambda] =1.5}, basic_4) , \ 

> [y(tau)] ,tau=-1.13. .1, [[y(0)=l]] , stepsize=0 . 01 , \ 

> axes=boxed,linecolor=red,arrows=none) :#expansion 

> p2 :=\ 

> DEplot (subs ({Omega [M] =0. 1, Omega [Lambda] =1.5}, \ 

> lhs(basic_4)=-rhs(basic_4)) , [y(tau)] ,tau=-3. .-1.14,\ 

> [[y(~3)=2.3]] ,stepsize=0.01,axes=boxed,linecolor=red,\ 

> arrows=none) : #contraction 

> display (pi, p2,title =< bouncing universe'); 
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4.3.8 Our universe (?) 



The modern observations of microwave background suggest JTTJ: £Ik ~ 
(flat), ~ 0.7, Qm ~ 0.3. The corresponding behavior of scaling fact 
looks as: 

> DEplot (subs ({Omega [M] =0.3, Omega [Lambda] =0 . 7} , \ 

> basic_4) , [y(tau)] ,tau=-.96. . .2, [[y(0)=l]] ,\ 

> stepsize=0 . 01 , axes=boxed, linecolor=red,arrows=none , \ 

> title= ( Our universe? '); 



Our universe? 




-1 -0.5 0.5 1 1.5 

lau 
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We are to note, that our comments about appropriate (or not appropriate) 
character of the considered models are not rigid. The picture of topologically 
nontrivial space-time can combine the universes with the very different local 
geometries and dynamics in the framework of a hypothetical Big Universe 
(" Tree of Universes" ) . 

4.4 Beginning 

4.4.1 Bianchi models and Mixmaster universe 

A good agreement of the isotropic homogeneous models with observations 
does not provide for their validity nearby singularity. In particular, in frame- 
work of analytical approach, we can refuse the isotropy of universe in the 
beginning of expansion (for review see jl4| ) . Let us consider the anisotropic 
homogeneous flat universe without rotation: 

> coord := [t, x, y, z] : 

> g_compts :=\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) :# metric components 

> g_compts [1 , 1] := -1: 

> g_compts [2,2] := a(t)"2: 

> g_compts [3,3] := b(t)"2: 

> g_compts [4,4] := c(t)~2: 

> g := create ( [-1,-1] , eval(g_compts) ) ; 

> ginv := invert ( g, 'detg' ): 



g := table( [index -char = 



[— 1, —1], compts 



-10 

a(t) 2 

h(t) 2 

c(f) 2 



]) 



The Einstein tensor is: 



100 



> # intermediate values 

> Dig := dlmetric( g, coord ) : 

> D2g := d2metric( Dig, coord ): 

> Cfl := Christoffell ( Dig ) : 

> RMN := Riemann( ginv, D2g, Cf 1 ) : 

> RICCI := Ricci( ginv, RMN ): 

> RS := Ricciscalar( ginv, RICCI ): 

> Estn := Einstein( g, RICCI, RS ):# Einstein tensor 

> displayGR(Einstein,Estn) ;# Its nonzero components 

The Einstein Tensor 



non — zero components : 

Gil = 

(m b(*)) (| ag» eg) + (| c(Q) (| a(t)) b(t) + (| eg)) (| b(t)) ag) 

a(t)b(t)c(t) 

= 

a(t) 2 gj| b(t)) eg) + c(Q) bg) + (| eg)) (| b(t))) 

b(t)c(i) 

G55 = 

b(t) 2 ((^ a(Q) eg) + (j| eg)) ag) + (| eg)) (| a(t))) 

a(t)c(t) 

= 

c(t) 2 gg ag)) bg) + (g bg)) ag) + (| b(t)) (| ag))) 

a(t)b(t) 

character : [—1,-1] 

In the vacuum state the right-hand sides of Einstein equations are 0. Then 

> E_eqn := get_compts (Estn) : 

> el := numer( E_eqn[l,l] ) = 0; 

> e2 :=\ 

> expand ( numer( E_eqn[2,2] )/a(t)~2 ) = 0; 

> e3 :=\ 

> expand ( numer( E_eqn[3,3] )/b(t)~2 ) = 0; 

> e4 :=\ 

> expand( numer( E_eqn[4,4] )/c(t)~2 ) = 0; 
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el := -(I b(t)) (| a(t)) c(f) - (f c(t)) (| a(t)) b(t)- 
(|c(i))(f b(t))a(t) = 



e£ := ( b(t)) c(t) + ( c(t)) b(t) + (| c(t)) (| b(t)) = 



e3 := (|1 a(t)) c(t ) + (|1 c(t)) a(t) + ( c(t) ) a(t) ) = o 



^ == a(t)) b(t) + b(t)) a(t) + (| b(t)) (| a(t)) = 



We will search the power-behaved solutions of this system ( Kasner metric ): 

> e5 := simplify (\ 

> subs({a(t)=a_0*t~p_l,b(t)=b_0*t~p_2,c(t)=c_0*t~p_3},\ 

> el) ); 

> e6 := simplify (\ 

> subs({a(t)=a_0*t~p_l,b(t)=b_0*t~p_2,c(t)=c_0*t~p_3},\ 

> e2) ); 

> e7 := simplify(\ 

> subs({a(t)=a_0*t~p_l,b(t)=b_0*t~p_2,c(t)=c_0*t~p_3},\ 

> e3) ); 

> e8 := simplify (\ 

> subs({a(t)=a_0*t~p_l,b(t)=b_0*t~p_2,c(t)=c_0*t~p_3},\ 

> e4) ); 

e5 := -b_0t(P- 2 - 2 +P- 1+ P-^ p.2 a_0 p.l c_0 

- c-0 t(P- 2 ~ 2 +P- 1 +P- 3 1 p_3 a_0 p.l b_0 

- c_0 tiP-Z-i+P-t+P- 3 ) p_3 b_0 p_2 a_0 = 

e6 := b_0t^- 3 - 2+ P-^ p_2 2 c_0 - b_0 t^- 3 - 2+ P-^ p_2 c_0 
+ c.O t (P- 3 - 2 +P- 2 ) p_3 2 b_0 - c_0 t (P- 3 ~ 2 +P- 2 ) p_3 b_0 
+ c.O t(P- 3 ~ 2 +P- 2 ) p_3 b.0 p-2 = 

el := a_0 j(p-S-2+j>-0 p _p c _ - a _0 £(p-s-2+p-0 pl c _ 
+ c_0 t { P- 3 - 2+ P-^ p_3 2 a_0 - c-0 t(P- 3 - 2 +P-^ p_3 a.O 
+ c.O *(p-s-2+p--0 p_3 a_0 p.l = 
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e8 := a_0t(P- 1 - 2+p - 2 '> p.l 2 b.O - a.O t^- 1 ' 2 ^-^ p.l b.O 
+ b.O -2+P-2) p _ 2 2 a _ Q _ b _ Q t (p_i-2+ P _2) p _ 2 a _ 

+ b.O jfr-'-a+P-*) p_2 a.O pA = 
If t is not equal to zero, we have: 

> e9 := factor( e5/t~ (p_2-2+p_l+p_3) ); 

> elO := factor ( e6/t" (p_3-2+p_2) ); 

> ell := factor ( e7/t~ (p_3-2+p_l) ); 

> el2 := factor ( e8/t" (p_2-2+p_l) ); 



e9 := - b_0 a_0 c.0 (p.2 pA + p.3 pA + p.3 p.2) = 



elO := b_0 c_0 {p.2 2 - p.2 + p.3 2 - p.3 + pS p.2) = 



ell := a_0 c_0 (pA 2 - pA + p.3 2 - p_3 + p_3 pA) = 



el2 := a_0 b.O (pA 2 - pA + p.2 2 - p.2 + p.2 pA ) = 

> el3 := e9/(-b_0*a_0*c_0) ; 

> el4 := elO/(b_0*c_0) ; 

> el5 := ell/(a_0*c_0) ; 

> el6 := el2/(a_0*b_0) ; 



el3 := p.2 p.l + p.3 p.l + p.3 p.2 = 



el4 := p.2 2 - p.2 + p.3 2 - p.3 + p.3 p.2 = 



el5 := p.l 2 -p.l + p.3 2 -p.3 + p.3 p.l =0 



el6 := p.l 2 - p.l + p.2 2 - p.2 + p.2 p.l = 



The following manipulation gives a connection between parameters: 
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> ((el4 + el5 + el6) - el3)/2; 

P-2 2 - p.2 + P-3 2 - pS + pA 2 - pA = 



That is pA 2 + p.2 2 + pS 2 = pA + p.2 + pS (3). 

> collect(el5-el6,{p_2,p_3}) ; 

> collect(el4-el5,{P-l>P-2}) ; 



p-3 2 + (pA - 1) p-3 - p.2 2 + (1 - pA )p.2 = 
P-2 2 + (p_3 - 1) p.2 - pA 2 + (1 - pS) pA = 

Hence 

pA (pA + - 1) = p.2 (p.2 + pA - 1), 
(pA + - 1) = pA (p-3 + p_i - 1). 
The last expressions suggest the simple solution: 

(p-3 + pA - 1) = p-2, (p-2 + pA - 1) = p-3 and 
+ p-2 - 1) = pA, (p-3 + - 1) = p-2, 
or 

+ pA - 1) = - (>_2 + - 1) = - and 
+ p_jg - 1) = - p_J, + pA - 1) = - 

The first results in 

> p_3 + p_l - 1 + p_2 + p_l - 1 + p_3 + p_2\ 

> - 1 - (p_l+p_2+p_3) ; 

p-3 + pA - 3 + p-2 

The second produces 
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> (p_3 + p_l - 1 + p_2 + p_l - 1 +\ 

> p_3 + p_2 - 1 + (p_l+p_2+p_3))/3; 

P-3 + pA - 1 + p.2 

These expressions in the combination with Eq.(3) suggest that there is only 
one independent parameter p: 

> solve({p_l + p_2 + p =\ 

> 1, p_l~2 + p_2~2 + p~2 = l},{p_l,p_2}): 

> sol_l := allvalues (°/ ) ; #f irst choice 

> solve({p_l + p_2 + p = 3, p_l~2 + p_2~2 + p~2 = 3},\ 
>{p_l,p_2}): 

> sol_2 := allvalues (°/ ) ; #second choice 

soU :={p_l=- l - P + l -- l -%l, p_2 = - l -p+ l -+ l -%l}, 

{pA = --p + - + - %1, p.2 = --p + - - - %1} 
%1 := ^-3p 2 + 2p + l 
soL2 := 

{p-2 = ~P+l + ^I{p-l)y/3, P-l =~P+|-^(P-1)V3}, 

So, if the parameters pA, p_2, pS are real, we have pA 2 + p_2 2 + pS 2 = 
pA + p-2 + p-3 =1 and 

> sol_p_l := factor(subs(sol_l,p_l)) ; 

> sol_p_2 := factor (subs (sol_l ,p_2) ) ; 

> sol_p_3 := 1 - H - %; 

soLpA := -Ip+i- I ^-(3^+1)^-1) 
sol.p.2 :=-\p+\ + \yJ-^P + l){jp-l) 



sol_p_3 := p 
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The values of the analyzed parameters are shown in the next figure: 



> plot({sol_p_l,sol_p_2,sol_p_3},p=-l/3. . 1,\ 

> axes=boxed, title= 'powers of t c ); 

powers of t 




-0.2 0.2 p0.4 0.6 0.8 



Hence, there are two directions of expansion and one direction of contraction 
for t-> oo (so-called "pancake" singularity) or one direction of expansion 
and two direction of contraction for t->0 ("sigar" singularity). 

Now let's consider more complicated situation with homogeneous anisotropic 
metric. 

We will use the tetrad notation of Einstein equations that allows to avoid 
the explicit manipulations with coordinates. In the synchronous frame we 
have for the homogeneous geometry (here Greek indexes change from 1 to 
3) Hi: 

90,0 = "I, 9a, /3 = Va,b {^a {e b )fj , 

where rj a ^ & is the time dependent matrix, {e a ) a is the frame vector (a changes 
from 1 to 3, that is the number of space triad). We will use the following 
representation for r] : 
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> eta :=\ 

> array (\ 

> [[a(t)-2,0,0] , [0,b(t)-2,0] , [0,0,c(t)-2]]);#a, b, c\ 

> are the time-dependent coefficients of anisotropic\ 

> deformation 

> eta_inv := inverse (eta) : 

> kappal : =\ 

> map(diff ,eta,t) : #kappal [a,b] =dif f (eta[a,b] ,t) 

> kappa2 := multiply(\ 

> map(dif f , eta,t) , eta_inv) ; #kappa2 [a] ~b=\ 

> dif f (eta[a,b] ,t)*eta_inv, i.e. we raise the index\ 

> by eta_inv 



a(t) 2 
b(t) 2 
c(t) 2 



k2 :-- 



a(t) 










b(t) 












c(t) 



The first vacuum Einstein equation has a form 

— . i iif " I k b K^ a — D 

— 2 dt a 4 a b ~ 



> eq_l :=\ 

> evalm( -trace (map(diff ,kappa2,t) ) /2 -\ 

> trace (multiply (kappa2 , transpose (kappa2) ) /4) )\ 

> = 0;#first Einstein equation 

^a(t) ^b(t) *cg) 

eg_i := , . qt , . . m , . — = 

* a(t) b(t) c(t) 



The next Einstein equations demand the definition of Bianchi structured 
constants: C ah = n ab + e abc a c , (C c ) ab = e abd C dc . Here e abc = e abc is 
the unit antisymmetric symbol, n ab is the symmetric "tensor", which can 
be expressed through its principal values iii , a c is the vector. The Bianchi 
types are: 
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Type 


a 


m 


n 2 


nz 


I 














II 





l 








VII 





l 


1 





VI 





l 


-1 





IX 





l 


1 


1 


VIII 





l 


1 


-1 


V 


1 











IV 


1 








1 


VII 


c 





1 


1 


in (c=i), vi (c + o) 


C 





1 


-1 



Table 1: Bianchi constants 



We will analyze type IX model, where C u = C 22 = C 33 = C23 1 = C31 
Ci 2 3 = 1. The curvature P a b = ±{2 C bd C ad + C db C ad + C bd C da - 

( {C b ) a + C a b )+ 5 a b [ (CV -2 C d f C df ]} is 



> Cab :=\ 

> array ([[1,0,0] , [0,1,0] , [0,0, 1] ] ) ; #CT (ab) 

> C_ab := multiply(eta,eta,Cab) ;#C[ab] 

> Ca_b := multiply(eta,Cab) ;#C[a] "b 

> C_a_b := multiply (Cab, eta) ;#C"b [a] ; 

> delta := array ([[1, 0,0] , [0,1,0] , [0,0,1]]) : 

> evalm( (4*multiply (Cab,C_ab) -\ 

> trace (C_a_b)*evalm(C_a_b+Ca_b) +\ 

> delta* (trace (C_a_b) ~2 -\ 

> 2*trace(multiply(Cab,C_ab))))/(2*det(eta))) :#here\ 

> we use the diagonal form of eta and symmetry\ 

> of structured constants 

> P := map(simplif y ,%) ; 



Cab : = 



1 
1 
1 



Cab :-- 



a(t) 4 
b(t) 4 
c(i) 4 
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Ca_b :-- 



a(t) 2 
b(t) 2 
c{t) 2 



CaJb :-- 



a(t) 2 
b(t) 2 
c(t) 2 



r i - a (t) 4 + b(t) 4 -2b(t) 2 c(t) 2 + c(t) 4 1 
2 a(t) 2 b(t) 2 c(t) 2 
p. n 1 -b(t) 4 + a(t) 4 -2a(t) 2 c(t) 2 + c(t) 4 
' 2 a(t) 2 b(t) 2 c(t) 2 

l c(t) 4 -a(t) 4 + 2a(t) 2 b(t) 2 -b(t) 4 
' 2 a(t) 2 b(t) 2 c(t) 2 

So, in the degenerate case a = b = c, t = const this model describes the 
closed spherical universe with a positive scalar curvature. 

Hence we can obtained the next group of vacuum Einstein equations: 

R " = -27m [(f vW) -Pa b ]=0 



> evalm(\ 

> -map(dif f , evalm(sqrt (det (eta) ) *kappa2) ,t)/\ 

> (2*sqrt (det (eta) ) ) ) : 

> evalm( map (simplify ,%) - P ); 
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(g aft)) bft) eg) + (| bft)) (| aft)) eft) + (| eft)) (| aft)) bft) 

aft) bft) eft) 

^(_ a(t) 4 + b(t) 4_ 2b(t) 2 c(t)2 + c(t) 4 ) 

+ 

aft) 2 bft) 2 eft) 2 ' ' 

(g, bft)) aft) eg) + (| bft)) (| aft)) eft) + (| eft)) (| bft)) aft) 

aft) bft) eft) 

(-bft) 4 + aft) 4 -2 aft) 2 eft) 2 + eft) 4 ) 
aft) 2 bft) 2 eft) 2 ' 



+ 



0, 0, 



(gr eg)) aft) bft) + (| eft)) (| aft)) bft) + (| eft)) (| bft)) aft) 



aft) bft) eft) 
1 eft) 4 - aft) 4 + 2 aft) 2 bft) 2 - bft) 4 



aft) 2 bft) 2 eft) 5 



The obtained expression gives three Einstein equations: 

> eq_2 :=\ 

> Diff( (diff (a(t) ,t)*b(t)*c(t)) ,t )/(a(t)*b(t)*c(t))\ 

> = ((b(t)-2-c(t)-2)~2 - a(t)"4 )/\ 

> (2*a(t)~2*b(t)~2*c(t)~2) ; 

> eq_3 :=\ 

> Diff( (diff (b(t) ,t)*a(t)*c(t)) ,t )/(a(t)*b(t)*c(t))\ 

> = ((a(t)-2-c(t)-2)~2 - b(t)~4 )/\ 

> (2*a(t)~2*b(t)~2*c(t)~2) ; 

> eq_4 :=\ 

> Diff( (diff (c(t) ,t)*b(t)*a(t)) ,t )/(a(t)*b(t)*c(t))\ 

> = ((a(t)-2-b(t)-2)~2 - c(t)"4 )/\ 

> (2*a(t)~2*b(t)~2*c(t)~2) ; 



(I a W) b W c ^ _ 1 (M*) 2 " eft) 2 ) 2 - aft) 4 



aft) bft) eft) 



aft) 2 bft) 2 eft) 5 



pn I (| bft)) aft) eft) _ 1 ( a ft) 2 - C ft) 2 ) 2 -bft) 4 



aft) bft) eft) 



aft) 2 bft) 2 cft) 5 
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_ |(|c(t))b(t)a(t) _ 1 (a(t) 2 -b(*) 2 ) 2 -c(t) 4 
69-4 " a(t)b(t)c(t) 2 a(t) 2 b(t) 2 c(t) 2 



The last group of Einstein equations results from 

Ra° = - \ K b C ( {C b ) ca - 5 a b (C d ) dc ) = 

It is obvious, that in our case this equation is identically zero (see expression 
for (C b ) ca ). The substitutions of a= e a , b= e 13 , c= e 7 (do not miss these 
exponents with frame vectors!), dt=a b c (d r ) and the expressions a = 

dr dt abc be dt 2 abc dr be abc i be \ be 

) a ] result in the simplification of the equations form: 



> eq_l :=\ 

> diff ((alpha(tau)+beta(tau)+gamma(tau)) , tau$2)/2 =\ 

> dif f (alpha(tau) , tau) *dif f (beta(tau) ,tau) +\ 

> dif f (alpha(tau) ,tau) *dif f (gamma(tau) ,tau) +\ 

> diff (beta(tau) ,tau) *dif f (gamma(tau) ,tau) ; 

> eq_2 :=\ 

> 2*diff (alpha(tau) ,tau$2) = (b(t)~2-c(t)~2)~2-a(t) ~4; 

> eq_3 :=\ 

> 2*diff (beta(tau) ,tau$2) = (a(t) ~2-c(t) ~2) ~2-b(t) ~4; 

> eq_4 :=\ 

> 2*diff (gamma (tau) ,tau$2) = (a(t)~2-b(t)~2)~2-c(t) ~4; 



1 ' d 2 , 1 / d 2 , 1 / d 



2 



eq.l := - a(r)) + - /3(r)) + - 7 (r)) = 
(Sf «W) (& fa)) + ( & a W) (& 7(r)) + (& /3(t)) (£ 7 (r)) 



e ^ := 2 (J^ a(r)) = (b(t) 2 - c(i) 2 ) 2 - a(i) 4 



eq_3 := 2 (3(r)) = (a(i) 2 - c(i) 2 ) 2 - b(t) 4 



eq.4 ■= 2 (J^ 7 (r)) = (a(t) 2 - b(t) 2 ) 2 - c(t) 4 



Comparison of obtained equations with above considered Kasner's type 
equations suggests that these systems are identical if the right-hand sides 
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of eq_2, eq_3 and eq_4 are equal to zero. Therefore we will consider the 
system's evolution as dynamics of Kasner's solution perturbation. In this 
case r = ln(i) + const (dt=a b c d r = t^ 1+P2+P3 ^ d r , as result d r = f 
) and J^a = pi , (5 = P2 , 7 = P3 • We assume that the right-hand 
sides are close to zero, but contribution of e^ a ^ - terms dominates. Hence: 

> eq_5 := diff (alpha(tau) , (tau,2)) =\ 

> -exp(4*alpha(tau) )/2; 

> eq_6 := diff (beta(tau) , (tau,2)) =\ 

> exp(4*alpha(tau) )/2; 

> eq_7 := diff (gamma(tau) , (tau,2)) =\ 

> exp(4*alpha(tau) )/2; 

eq.5 :=^- 2 *(T) = - 1 -e^)) 



eq.6 : = A_/3( r ) = I e (^W) 



The first equation is analog of equation describing the one-dimensional mo- 
tion ( a is the coordinate) in the presence of exponential barrier: 

> p:=\ 

> dsolve({eq_5,alpha(5)=-2,D(alpha) (5)=-0.5},alpha(tau) ,\ 

> type=numeric) : 

> with(plots) : 

> odeplot (p, [tau,dif f (alpha (tau) ,tau)] , -5 . . 5 , \ 

> color=green) : 

> plot(exp(4*(-0.5)*tau)/2,tau=-5. . 5 , color=red) : 

> display (°/ , °/,°h, view=- 1 . . l,axes=boxed,\ 

> title =< reflection on barrier'); 
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We can see, that the" particle" with initial velocity p\ = -0.5 (green curve) 
is reflected on barrier (red curve) with change of velocity sign. But eg_5, 
eq_6, eq_7 results in p\ + P2 = const, p\ + ps = const therefore 

(p«««) a = p 2 + Pl - (p^)! = p 2 + 2 Pl , (p new ) 3 = P3 + Pl - (;p new )i = 

P3 + 2 pi , 

t = ab c = e« 1 + 2 P 1 ) T ) , 

and the result is: 

i vi \ ( P2+ 2 Pi \ (E2+lZL\ 

a = t ( 1 + i Pi J , b = V 1+2 pi ' , c = V 1+2 "i ' 

The following procedure gives a numerical algorithm for iteration of powers 
of t. 
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> iterations := proc(iter ,p) 

> p_l_old :=\ 

> evalhf ( -l/2*p+l/2-l/2* S qrt(-(3*p+l)*(p-D) ); 

> p_2_old :=\ 

> evalhf ( -l/2*p+l/2+l/2*sqrt(-(3*p+l)*(p-l)) ); 

> p_3_old := evalhf ( p ): 

> for m from 1 to iter do 

> if (p_l_old<0) then 

> pp := evalhf (l+2*p_l_old) : 

> p_l_n := evalhf (-p_l_old/pp) : 

> p_2_n := evalhf ( (p_2_old+2*p_l_old)/pp) : 

> p_3_n := evalhf ( (p_3_old+2*p_l_old) /pp) : 

> else fi: 

> if (p_2_old<0) then 

> pp := evalhf (l+2*p_2_old) : 

> p_l_n := evalhf ( (p_l_old+2*p_2_old) /pp) : 

> p_2_n := evalhf (-p_2_old/pp) : 

> p_3_n := evalhf ( (p_3_old+2*p_2_old) /pp) : 

> else fi: 

> if (p_3_old<0) then 

> pp := evalhf (l+2*p_3_old) : 

> p_l_n := evalhf ((p_l_old+2*p_3_old)/pp) : 

> p_2_n := evalhf ( (p_2_old+2*p_3_old) /pp) : 

> p_3_n := evalhf (-p_3_old/pp) : 

> else fi: 

> p_l_old := p_l_n: 

> p_2_old := p_2_n: 

> p_3_old := p_3_n: 

> if m = iter then pts :=\ 

> [p_l_old, p_2_old,p_3_old] fi; 

> od: 

> pts 

> end: 

> with(plots) : 

> pointplot3d(\ 

>{seq(iterations(i,0.6) , i=l . . 8) } , symbol=diamond, \ 

> axes=B0XED,labels=[p_l,p_2,p_3] ) ; 

> pointplot3d(\ 

>{seq(iterations(i,0.6) , i=7 . . 12)},symbol=diamond,\ 

> axes=B0XED,labels=[p_l,p_2,p_3] ) ; 



114 




One can see, that the evolution has character of switching between different 
Kasner's epochs (all points lie on the section of sphere pA 2 + p_2 2 + pS 2 
= 1 by plane pA + p_2 + pS =1). The negative power of t increases 
the corresponding scaling factor if f->0. In the upper figure the periodical 
change of signs takes place between y- and z- directions. The universe 
oscillates in these directions and monotone contracts in x-direction (f->0). 
Next the evolution changes (lower figure): x-z oscillations are accompanied 
by the contraction in y-direction. The whole picture is the switching between 



115 



different oscillations as the result of logarithmical approach to singularity. It 
is important, that there is the strong dependence on initial conditions that 
can cause the chaotical scenario of oscillations: 



> pointplot3d(\ 

>{seq( iterations (i ,0 . 601) , i=l . . 100) } , symbol=diamond, \ 

> axes=B0XED,labels=[p_l,p_2,p_3] ,\ 

> title= c Kasner ' s epochs (I)'); 

> pointplot3d(\ 

>{seq( iterations (i ,0 . 605) , i=l . . 100) } , symbol=diamond, \ 

> axes=B0XED,labels=[p_l,p_2,p_3] A 

> title= c Kasner ' s epochs (II)'); 
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So, we can conclude that the deflection from isotropy changes the scenario 
of the universe's evolution in the vicinity of singularity, so that the dynamics 
of the geometry looks like chaotical behavior of the deterministic nonlinear 
systems (so-called Mixmaster universe , see [jl6|). 



4.4.2 Inflation 

The above-considered standard homogeneous models face the challenge of 
so-called horizon . The homogeneity of the universe results from the causal 
connection between its different regions. But the distance between these 
regions is defined by the time of expansion. Let's find the maximal distance 
of the light signal propagation in the expanding region. As ds 2 = the 
accompanying radial coordinate of horizon is 

> r = Int(l/a(t) ,t=0. .tO) ;# tO is the age\ 

> of the universe 

> #or 

> r = Int(l/y(tau) ,tau=-tO*HO. .0)/(aO*HO) ; 

to i 

dt 



a(t) 



to hoy{t) 



aO HO 



For instance, in the Einstein-de Sitter model we have 

> int(\ 

> subs(y(tau)=l/4*(12*tau+8)~(2/3) ,l/y(tau)) ,\ 

> tau=-2/3 . . 0) / (aO*HO) : # we use the above obtained\ 

> expressions for this model 

> simplify( subs (a0=l/4* (12*0+8) ~ (2/3) ,%) .radical ); 

2 w 



and in the Friedmann-Lemaitre: 
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> int(\ 

> subs (y (phi) =Omega [M] * (1-cos (2*phi) ) /\ 

> (2*b) ,l/2*0mega[M]*(2-2*cos(2*phi))/\ 

> (b"(3/2))/y(phi)) ,phi )/aO/HO;# we use \ 

> the above obtained expressions for this model 

> simplify( subs (K=-l , subs (b=-K/aO~2/H(T2 //„)), \ 

> radical , symbolic) ;#b=0mega[M] -l=0mega[K] and K=-l\ 

> (spherical universe) 

VbaOHO 



2(f) 

The last example is of interest. The increase of 4> from to ir /2 corresponds 
to expansion of universe up to its maximal radius (see above). At the mo- 
ment of maximal expansion the observer can see the photons from antipole 
of universe that corresponds to formation of full causal connection: 

> plot( [2*phi-sin(2*phi) ,2*phi] ,\ 

> phi=0..Pi, title='age of universe and time of\ 

> photon travel (a.u.)'); 

age of universe and time of photon travel (a.u.) 




05 1.5. IT 2.5 

phi 



This figure shows the age of Friedmann-Lemaitre universe in comparison 
with the time of photon propagation from outermost point of universe. For 
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(ft < it /2 there are the regions, which don't connect with us. After this 
points there is the causal connection. At the moment of recollapse the pho- 
tons complete revolutions around universe. 



The different approaches were developed in order to overcome the problem 
of horizon (see above considered loitering and Mixmaster models). But the 
most popular models are based on the so-called inflation scenario , which can 
explain also the global flatness of universe (see, for example, [17], on-line re- 
views HI, HI). 



For escape the horizon problem we have to suppose the accelerated expan- 
sion with Jp- a(i) > 0. As result of this condition, the originally causally- 
connected regions will expand faster, than the horizon expands. So, we ob- 
serve the light from the regions, which was connected at the early stage of 
universe's evolution. In the subsection "Standard models" we derived from 
Einstein equations the energy conservation low basic_2: 

4^ = -4vr(p(t ) + 4)) + A 

This equation demonstrates that in the usual conditions {p > 0) the pressure 
is a source of gravitation, but if p < - | (A = 0)orA>4-7r/9 (p=0) the 
repulsion dominates and the expansion is accelerated. As the source of the 
repulsion, the scalar "inflaton" field (ft is considered. 



Let us introduce the simple Lagrangian for this homogeneous field: 

L-lirJ^tftg* -V(<j>), 

where V is the potential energy. The procedure, which was considered in 
first section, allows to derive from this Lagrangian the field equation: 




In the case of flat RW metric this equation results in: 
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> coord := [t, r, theta, phi]: 

> g_compts : =\ 

> array (symmetric , sparse , 1 . .4, 1 . .4) : 

> g_compts [1 , 1] := -1: 

> g_compts [2,2] := a(t)~2: 

> g_compts [3,3] := a(t)~2*r~2: 

> g_compts [4,4] := a(t) ~2*r~2*sin(theta) ~2 : 

> g := create ( [-1,-1] , eval (g_compts) ) : #def inition\ 

> of flat RW metric 

> dlg: = dlmetric(g, coord): 

> g_inverse:= invert ( g, ' detg'): 

> Cfl:= Christoffell ( dig ) : 

> Cf2:= Christoffel2 ( g_inverse, Cfl ): 

> det_scal := simplify (\ 

> sqrt( -det( get_compts (g) )) .radical, symbolic ): 

> g_det_sq :=\ 

> create ([], det_scal) : #sqrt (-g) 

> field :=\ 

> create ([], phi (t) ) : #f ield 

> #calculation of first term in the field equation 

> cov_diff( field, coord, Cf2 ): 

> T := prod(g_inverse,%) : 

> contract (T, [2,3] ) : 

> Tt := prod(g_det_sq,°/„) : 

> Ttt := cov_diff( Tt, coord, Cf2 ): 

> get_compts (%) [1 , 1] /det_scal : 

> expand (-°/ ) : 

> field_eq := 7, + dif f (V(phi) ,phi) = 0; 

fields := 3 — + {— 2 m + (- V(0) = 



Now let's define the energy-momentum tensor for first Einstein equation: 

Ty = gij L - gji f x _i -qj—^ L , 

where f x _ { = ■£:(/>■ As consequence, p = [- \ <j)) 2 + V( (f> )] + (f)) 2 = 
\ (I +V(h, and (for i, j = 1,2,3) T y = 9ij ( \ (| 4>f - V( <j> )) => p 
= i(^0) 2 - V( <p). Last expression emplies that the necessary condition 
p < - | may result from large V. 

The dynamics of our system is guided by coupled system of equations: 

> basic_inf_l :=\ 

> K/(a(t)"2)+diff (a(t) ,t)~2/(a(t)~2) =\ 

> 8/3*Pi*(l/2*diff (phi(t) ,t)~2+V(phi)) ;#see above E[l] 

> basic_inf_2 := field_eq; 
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fa S!CJn /.2 := 3 lAM»Aig» + ( J| + (i V («) = 



The meaning of K is considered earlier, here we will suppose K=0 (local 
flatness). In the slow-rolling approximation J^- <f>(t) , ( <fi(t)) 2 ->0 and for 
potential V( (f> ) = ^j- we have: 



> basic_inf _l_n :=\ 

> subs(K=0,K/(a(t)~2)+diff (a(t) ,t)~2/(a(t)"2)) =\ 

> 8/3*Pi*m~2*phi(t)~2/2;#see above E[l] 

> basic_inf _2_n := op(l , lhs (f ield_eq) ) +\ 

> subs (phi=phi (t) ,expand(subs(V(phi)=m~2*phi~2/2,\ 

> op(3,lhs(field_eq)))))=0; 



basic Jnf.l.n := idt ^ = -Trm 2 ^) 2 
a(t) 2 3 v ; 



basic.mf-2-n := 3 a(t)) + m 2 <f>(t) = 

a(t) 

> dsolve({basic_inf _l_n,basic_inf _2_n, \ 

> phi(0)=phiO,a(0)=aO},{phi(t) ,a(t)}) : 

> expand (%) ; 



I W e - 2/34-0 y^mV3t ' ^ W 6 y/lT I 



One can see, that for <fio >>1 there is an quasi-exponential expansion (in- 
flation) of universe: 



a(t)= a e< Ht ) , where H = 2 ^ <p 



m 



When t ~ 2 ^- fa the regime of slow rolling ends and the potential energy 
of inflaton field converts into kinetic form that causes the avalanche-like 
creation of other particles from vacuum (so-called "reheating"). Only at 
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this moment the scenario become similar to Big Bang picture. What is the 
size of universe after inflation? 

> a(t)=subs(\ 

>{H=2*sqrt (Pi/3) *phi [0] *m , t=2*sqrt (3*Pi) /m*phi [0] } , \ 

> a[0]*exp(H*t)) ; 

a(t) = a e (4 "^ 2) 

If at the beginning of inflation the energy and size of universe are defined by 
Plank density and length, then m 2 (0 2 )o "1, «o ~ 10^~ 33 ^ cm. The necessary 
value of m ~ ' 10^~^ is constrained from the value of fluctuation of back- 
ground radiation ( ~10^ 5) ). As result a ~ lO^" 33 ) e ( ^ ) = lO^ 33 ) e ( 47rl ° 12 ) 
cm , that is much larger of the observable universe. This is the reason why 
observable part of universe is flat, homogeneous and isotropic. 

Let's return to basic_inf-2. The reheating results from oscillation of inflaton 
field around potential minimum that creates new boson particles. The last 
damps the oscillations that can be taking into consideration by introducing 
of dumping term in basic-inf_2: 

> reh_inf :=\ 

> 3*diff (phi(t) ,t)*H+Gamma*diff (phi(t) ,t) +\ 

> diff (phi(t) , (t,2)) = -m"2*phi(t) ; 

reh.inf := 3 (|- <£(*)) H + + = 



Here H is the Hubble constant, T is the decay rate, which is defined by 
coupling with material boson field. When H > T , there are the coherent 
oscillations: 

> dsolve(subs(Gamma=0,reh_inf ) ,phi(t)) ; 

> diff(°/.,t); 

4>(t) = _C1 e (- 1/2 (3 H+y/(3 H-2 m) (3 H+2 m))t) + 
C2 e^ 1 / 2 (3#-\/( 3 - H_2m ) (3H+2m)) t) 
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§- t cf>(t)=-Cl {~\h- I%l) e (- 1 /2(3i?+%l)t) + 

_^(-^+I%l) e (-l/2(3H-%l) < ) 

%1 := ^(3#-2m)(3F + 2m) 

The last expression demonstrates the decrease of energy density via time 
increase ( p ~ (-§i4>) 2 ) and, as result, the decrease of H (see basic-infA) . 
This leads to H < T , that denotes the reheating beginning. 

In the conclusion of this short and elementary description it should be noted, 
that the modern inflation scenarios modify the standard model essentially. 
The Universe here looks like fractal tree of bubbles-universes, inflating and 
recollapsing foams, topological defects etc. The main features of this picture 
wait for thorough investigations. 



5 Conclusion 



Our consideration was brief and I omitted some important topics like per- 
turbation theory for black holes and cosmological models, Penrose diagrams 
and conformal mappings for investigation of causal structure, gravitational 
waves etc. But I suppose to advance this worksheet hereafter. 
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